Regression Project

In this project, your goal is to build regression models of housing prices. The models should learn from data and be able to predict the median house price in a district (which is a population of 600 to 3000 people), given some predictor variable

Das Vorgehen in diesem Notebook orientiert sich am Data Science Lifecycle. Es wird eine Modellpräzision von > 70% gewünscht. Das bedeutet, dass der R-squared Wert des besten Modells >0,7 sein soll.

Python Setup

Zunächst werden mit der Import Methode wichtige Python Bibliotheken importiert.

import pandas as pd
import numpy as np 

import seaborn as sns
sns.set_theme(style='ticks')

import matplotlib.pyplot as plt
%matplotlib inline

import statsmodels.formula.api as smf
import statsmodels.api as sm
import statsmodels.gam.api as smg

from sklearn.compose import ColumnTransformer
from sklearn.compose import make_column_selector as selector
from sklearn.pipeline import Pipeline
from sklearn.impute import SimpleImputer
from sklearn.preprocessing import StandardScaler, OneHotEncoder
from sklearn.model_selection import train_test_split
import sklearn.metrics as skm
from sklearn.linear_model import LinearRegression
from sklearn.linear_model import Lasso
from sklearn.preprocessing import SplineTransformer
from sklearn.linear_model import Ridge
from sklearn.pipeline import make_pipeline

import plotly.express as px

Datenaufbereitung

Daten laden

Importieren der Daten aus der CSV in Github und Erstellung eines Pandas Dataframes zur Bearbeitung und Verwendung der Daten in diesem Notebook.

GIT = "https://raw.githubusercontent.com/jan-kirenz/project-annikas428/main/"
DATA = "project_data.csv"
TOKEN = "?token=GHSAT0AAAAAABPURMBOEW2WS3JEBXGI5TPAYPJRBIQ"

df_initial = pd.read_csv(GIT + DATA + TOKEN)

Inspect Data

Im ersten Schritt ist es wichtig, sich einen Überblick über die Daten zu verschaffen. Begonnen wird mit einer Anzeige einiger Datensätze.

df_initial
longitude latitude housing_median_age total_rooms total_bedrooms population households median_income median_house_value ocean_proximity price_category
0 -122.23 37.88 41.0years 880 129.0 322 126 8.3252 452600.0$ NEAR BAY above
1 -122.22 37.86 21.0 7099 1106.0 2401 1138 8.3014 358500.0 NEAR BAY above
2 -122.24 37.85 52.0 1467 190.0 496 177 7.2574 352100.0 NEAR BAY above
3 -122.25 37.85 52.0 1274 235.0 558 219 5.6431 341300.0 NEAR BAY above
4 -122.25 37.85 52.0 1627 280.0 565 259 3.8462 342200.0 NEAR BAY above
... ... ... ... ... ... ... ... ... ... ... ...
20635 -121.09 39.48 25.0 1665 374.0 845 330 1.5603 78100.0 INLAND above
20636 -121.21 39.49 18.0 697 150.0 356 114 2.5568 77100.0 INLAND above
20637 -121.22 39.43 17.0 2254 485.0 1007 433 1.7000 92300.0 INLAND above
20638 -121.32 39.43 18.0 1860 409.0 741 349 1.8672 84700.0 INLAND above
20639 -121.24 39.37 16.0 2785 616.0 1387 530 2.3886 89400.0 INLAND above

20640 rows × 11 columns

Beim Sichten der Datensätze fällt auf, dass im ersten Datensatz zwei Unstimmigkeiten im Vergleich zu den anderen Datensätzen vorliegen:
median_house_value ist nur hier mit der Währung Dollar angegeben und hinter housing_median_age steht der Zusatz “years”.

Nun folgt Pandas Dataframe Info Methode. Diese zeigt die Variablen und deren Datentypen, sowie die Anzahl der non-null values je Spalte an.

df_initial.info()
<class 'pandas.core.frame.DataFrame'>
RangeIndex: 20640 entries, 0 to 20639
Data columns (total 11 columns):
 #   Column              Non-Null Count  Dtype  
---  ------              --------------  -----  
 0   longitude           20640 non-null  float64
 1   latitude            20640 non-null  float64
 2   housing_median_age  20640 non-null  object 
 3   total_rooms         20640 non-null  int64  
 4   total_bedrooms      20433 non-null  float64
 5   population          20640 non-null  int64  
 6   households          20640 non-null  int64  
 7   median_income       20640 non-null  float64
 8   median_house_value  20640 non-null  object 
 9   ocean_proximity     20640 non-null  object 
 10  price_category      20640 non-null  object 
dtypes: float64(4), int64(3), object(4)
memory usage: 1.7+ MB

Die Variablen housing_median_age und und median_house_value wurden als Datentyp object definiert. Dies muss genauer untersucht werden, da bei einer Jahreszahl bzw. einem monetären Wert ein numerischer Datentyp (int oder float) angemessener ist. Die Ursache hierfür ist auf die Zusatzangaben “$” und “years” zurückzuführen.
Die Varialben ocean_proximity und price_category sind keine numerischen sondern kategoriale Datentypen. Diese lassen sich jeweils in einer internen Rangordnung anordnen. Sie sind daher ordinal. Der richtige Datentyp hierfür ist category. Die Transformation der Datentypen erfolgt in Kapitel “Data Transformation”.
Neben den Datentypen wird auch die Anzahl der non-null values je Variable angezeigt. Hier haben alle Variablen eine Anzahl von 20640 non-null values. Lediglich die Variable total_bedrooms liegt unter diesem Wert. Dies lässt darauf schließen, dass hier einige Datensätze keinen Wert haben. Dies gilt es im Folgenden genauer zu untersuchen.

# show missing values -if present - in yellow
sns.heatmap(df_initial.isnull(), 
            yticklabels=False,
            cbar=False, 
            cmap='viridis');

print(df_initial.isnull().sum())
longitude               0
latitude                0
housing_median_age      0
total_rooms             0
total_bedrooms        207
population              0
households              0
median_income           0
median_house_value      0
ocean_proximity         0
price_category          0
dtype: int64
_images/regression_13_1.png

Wie bereits vermutet, enthält die Variable total_bedrooms 207 null values. Die Abbildung visualisiert diese nicht befüllten Datensätze und deren Verteilung im Datenset in gelb. Die null-values sind zufällig im Datenset verteilt. Da das Datenset über 20000 Datensätze enthält, werden die Zeilen mit den null-values im Kapitel Transform Data entfernt. Alternativ wäre auch eine Belegung der null-values beispielsweise mit dem Durchschnittswert der Variable möglich. Darauf wird hier jedoch verzichtet, genug andere Datensätze vorliegen.

Datentransformation

Entfernung der null values

# drop rows with missing values
df = df_initial.dropna()
print(df.isnull().sum())
longitude             0
latitude              0
housing_median_age    0
total_rooms           0
total_bedrooms        0
population            0
households            0
median_income         0
median_house_value    0
ocean_proximity       0
price_category        0
dtype: int64

Entfernung unwichtiger Features

Das Feature -price_category_ wird entfernt, da es sich direkt aus der Variablen median_house_value ableitet und daher nicht zur Berechnung des Wertes genutzt werden darf. Dieses Feature wird im Notebook Classification als abhängige Variable genutzt.

#drop irrelevant features
df = df.drop(['price_category'], axis=1)

Konvertierung von Datentypen

Die mithilfe der info() Methode identifizierten Features, welche eine Datentyp Konvertierung benötigen, werden nun bearbeitet.

# Convert data types Ocean_Proximity and Price_Category to Data Type "Category" as they are categorical variables (ordinal)
df.ocean_proximity = df.ocean_proximity.astype('category')

# Housing_Median_Age and Median_House_Value are actually int64 values, but pandas did not get this because of signs like '€'
df.housing_median_age = df.housing_median_age.str.replace('years', '', regex=False)
df.median_house_value = df.median_house_value.str.replace('$', '', regex=False)
df.housing_median_age = pd.to_numeric(df.housing_median_age, errors='coerce').astype('int64')
df.median_house_value = pd.to_numeric(df.median_house_value, errors='coerce').astype('float64')
df.total_bedrooms = df.total_bedrooms.astype('int64')
df.info()
df.head()
<class 'pandas.core.frame.DataFrame'>
Int64Index: 20433 entries, 0 to 20639
Data columns (total 10 columns):
 #   Column              Non-Null Count  Dtype   
---  ------              --------------  -----   
 0   longitude           20433 non-null  float64 
 1   latitude            20433 non-null  float64 
 2   housing_median_age  20433 non-null  int64   
 3   total_rooms         20433 non-null  int64   
 4   total_bedrooms      20433 non-null  int64   
 5   population          20433 non-null  int64   
 6   households          20433 non-null  int64   
 7   median_income       20433 non-null  float64 
 8   median_house_value  20433 non-null  float64 
 9   ocean_proximity     20433 non-null  category
dtypes: category(1), float64(4), int64(5)
memory usage: 1.6 MB
longitude latitude housing_median_age total_rooms total_bedrooms population households median_income median_house_value ocean_proximity
0 -122.23 37.88 41 880 129 322 126 8.3252 452600.0 NEAR BAY
1 -122.22 37.86 21 7099 1106 2401 1138 8.3014 358500.0 NEAR BAY
2 -122.24 37.85 52 1467 190 496 177 7.2574 352100.0 NEAR BAY
3 -122.25 37.85 52 1274 235 558 219 5.6431 341300.0 NEAR BAY
4 -122.25 37.85 52 1627 280 565 259 3.8462 342200.0 NEAR BAY

Mithilfe der string replace Methode werrden die zuvor identifizierten Zeichen “years” und “$” entfernt. Anschließend ist eine Transformation der Features in integer bzw. float Werte möglich. Für housing_median_age wird integer gewählt, da die Jahreszahlen als Ganze Zahlen also discrete Werte angegeben sind.
Für median_house_value wird der Datentyp float benötigt, da die Ergebnisse einer linearen Regression stehts continous sind.

Hinzufügen neuer Features

df = df.assign(people_per_household=lambda df: df.population / df.households)
df = df.assign(rooms_per_household=lambda df: df.total_rooms / df.households)
df = df.assign(bedrooms_per_people=lambda df: df.total_bedrooms / df.population)
df = df.assign(bedrooms_per_rooms=lambda df: df.total_bedrooms / df.total_rooms)

Diese Features werden dem Datenset hinzugefügt, in der Hoffnung, dass diese für die Verwendung in den Modellen nützlich sein könnten. Untersucht wird dies im Kapitel EDA.

Datensplit

Split des Datensets in 80% Trainingsdaten und 20% Testdaten. Die Testdaten werden anschließend bis zur finalen Evaluation der besten Modelle nicht verwendet.

train_dataset = df.sample(frac=0.8, random_state=0) 
test_dataset = df.drop(train_dataset.index)

train_dataset
longitude latitude housing_median_age total_rooms total_bedrooms population households median_income median_house_value ocean_proximity people_per_household rooms_per_household bedrooms_per_people bedrooms_per_rooms
14185 -117.08 32.70 37 2176 418 1301 375 2.8750 98900.0 NEAR OCEAN 3.469333 5.802667 0.321291 0.192096
6125 -117.91 34.11 20 3158 684 2396 713 3.5250 153000.0 <1H OCEAN 3.360449 4.429173 0.285476 0.216593
14095 -117.10 32.75 11 2393 726 1905 711 1.3448 91300.0 NEAR OCEAN 2.679325 3.365682 0.381102 0.303385
14359 -117.22 32.74 52 1260 202 555 209 7.2758 345200.0 NEAR OCEAN 2.655502 6.028708 0.363964 0.160317
18004 -121.99 37.29 32 2930 481 1336 481 6.4631 344100.0 <1H OCEAN 2.777547 6.091476 0.360030 0.164164
... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
6168 -117.93 34.09 35 1891 353 1093 382 4.0167 165500.0 <1H OCEAN 2.861257 4.950262 0.322964 0.186674
4174 -118.21 34.10 36 2000 533 1234 535 3.7437 241700.0 <1H OCEAN 2.306542 3.738318 0.431929 0.266500
3593 -118.48 34.24 32 2621 412 1285 414 6.6537 267600.0 <1H OCEAN 3.103865 6.330918 0.320623 0.157192
5366 -118.39 34.04 44 1873 286 635 283 5.5951 461300.0 <1H OCEAN 2.243816 6.618375 0.450394 0.152696
17879 -122.00 37.40 17 5121 1017 1470 968 2.9706 81300.0 <1H OCEAN 1.518595 5.290289 0.691837 0.198594

16346 rows × 14 columns

Deskriptive Statistik

Kennzahlenanalyse numerischer Variablen

Mithilfe der describe() Methode werden die wichtigsten Kennzahlen numerischer Variablen (bspw. Durchschnitt, min und max Werte) ausgegeben.

#summary statistics for all numerical columns
train_dataset.describe().T
count mean std min 25% 50% 75% max
longitude 16346.0 -119.564154 2.002618 -124.350000 -121.790000 -118.490000 -118.000000 -114.470000
latitude 16346.0 35.630318 2.138574 32.550000 33.930000 34.250000 37.710000 41.950000
housing_median_age 16346.0 28.664505 12.556764 1.000000 18.000000 29.000000 37.000000 52.000000
total_rooms 16346.0 2622.235776 2169.548287 11.000000 1448.000000 2119.000000 3120.750000 39320.000000
total_bedrooms 16346.0 535.281659 418.469078 3.000000 296.000000 432.500000 644.000000 6445.000000
population 16346.0 1416.087055 1103.842065 3.000000 784.250000 1164.000000 1711.000000 28566.000000
households 16346.0 496.758167 379.109535 3.000000 280.000000 408.000000 600.000000 6082.000000
median_income 16346.0 3.869337 1.902228 0.499900 2.555675 3.533200 4.744225 15.000100
median_house_value 16346.0 206916.154411 115676.394484 14999.000000 119300.000000 179700.000000 265900.000000 500001.000000
people_per_household 16346.0 3.089529 11.525259 0.692308 2.428850 2.816514 3.282584 1243.333333
rooms_per_household 16346.0 5.423747 2.160962 0.846154 4.443697 5.230149 6.051703 62.422222
bedrooms_per_people 16346.0 0.403370 0.223233 0.000670 0.315480 0.372646 0.443446 8.750000
bedrooms_per_rooms 16346.0 0.213188 0.057891 0.100000 0.175591 0.203263 0.239668 1.000000

Der Durchschnitt (mean) der Variablen housing_median_age über alle Datensätze beträgt 28,64 Jahre. Die Standardabweichung, also die durchschnittliche Abweichung der Werte vom Durchschnitt, beträgt 12,56 Jahre. Das Minimum liegt bei 1 Jahr, das Maximum der 52 Jahren. 25% der Datensätze liegen unter 18 Jahren (1. Quartil). 50% zwischen unter 29 Jahren (2. Quartil) und 75% zwischen unter 37 Jahren (3. Quartil).
Im Durchschnitt beträgt der Wert people_per_household 3,09. Die durchschnittliche Abweichung von diesem Wert beträgt 11,53. Das Minimum liegt bei 0,69 und das Maximum bei 1243 Personen pro Haushalt. 25% der Datensätze liegen zwischen 0,69 und 2,43. 50% unter 2,82 und 75% der Datensätze unter 3,28 Personen pro Haushalt.

print(train_dataset.quantile(0.99, interpolation='nearest'))
longitude                 -116.290000
latitude                    40.630000
housing_median_age          52.000000
total_rooms              11021.000000
total_bedrooms            2208.000000
population                5723.000000
households                1955.000000
median_income               10.557500
median_house_value      500001.000000
people_per_household         5.286920
rooms_per_household         10.500000
bedrooms_per_people          0.941606
bedrooms_per_rooms           0.407472
Name: 0.99, dtype: float64

Die Quantile() Methode gibt für jede Variable einen Wert an, unter welchem (in diesem Beispiel) 99% der Datensätze liegen.
Beispielsweise liegen 99% der Datensätze unter einer Anzahl total_rooms in Höhe von 11021.
99% der Distrikte haben eine population unter 5723 Einwohnern.
Im Nachfolgenden werden einige mögliche Ausreiser identifiziert.

Dieser Datensatz steht für ein Distrikt, welches aus nur 1 Haushalt besteht. In diesem Haushalt leben 13 Leute, die sich 8 Räume, wovon einer ein Schlafzimmer ist, teilen. Der Wert des Hauses beträgt 500001 $.
Dies erscheint etwas unwahrscheinlich. Womöglich handelt es sich um eine Falscheingabe. Der Datensatz wird daher entfernt.

train_dataset.loc[train_dataset.people_per_household > 100]
longitude latitude housing_median_age total_rooms total_bedrooms population households median_income median_house_value ocean_proximity people_per_household rooms_per_household bedrooms_per_people bedrooms_per_rooms
19006 -121.98 38.32 45 19 5 7460 6 10.2264 137500.0 INLAND 1243.333333 3.166667 0.000670 0.263158
16669 -120.70 35.32 46 118 17 6532 13 4.2639 350000.0 NEAR OCEAN 502.461538 9.076923 0.002603 0.144068
3364 -120.51 40.41 36 36 8 4198 7 5.5179 67500.0 INLAND 599.714286 5.142857 0.001906 0.222222

Mehr als 50 Leute in einem Haushalten wirken sehr unrealistisch. 99% der Datensätze liegen unter einem Wert von 5,39 Personen pro Haushalt. Daher werden diese Datensätze als Ausreiser identifiziert.

train_dataset.loc[train_dataset.bedrooms_per_people > 8]
longitude latitude housing_median_age total_rooms total_bedrooms population households median_income median_house_value ocean_proximity people_per_household rooms_per_household bedrooms_per_people bedrooms_per_rooms
11862 -121.25 40.27 25 958 245 28 16 2.625 67500.0 INLAND 1.75 59.875 8.75 0.255741

Auch diese Werte liegen weit über dem 99% Perzentil der Varialben bedrooms_per_people und werden daher als Ausreiser entfernt. Zwei dieser Werte sind auch bereits rooms_per_household Ausreiser.

train_dataset.loc[train_dataset.rooms_per_household > 55]
longitude latitude housing_median_age total_rooms total_bedrooms population households median_income median_house_value ocean_proximity people_per_household rooms_per_household bedrooms_per_people bedrooms_per_rooms
12447 -114.49 33.97 17 2809 635 83 45 1.6154 87500.0 INLAND 1.844444 62.422222 7.650602 0.226059
11862 -121.25 40.27 25 958 245 28 16 2.6250 67500.0 INLAND 1.750000 59.875000 8.750000 0.255741
1913 -120.06 39.01 19 2967 528 112 48 4.0714 437500.0 INLAND 2.333333 61.812500 4.714286 0.177958
#drop people per household outliers
train_dataset = train_dataset.drop([3364, 16669, 19006])
#drop bedroooms per people outliers
train_dataset = train_dataset.drop([11862])
#drop rooms per household_outliers
train_dataset = train_dataset.drop([1913, 12447])
#summary statistics for all numerical columns
train_dataset.describe().T
count mean std min 25% 50% 75% max
longitude 16340.0 -119.564056 2.002423 -124.350000 -121.790000 -118.490000 -118.000000 -114.470000
latitude 16340.0 35.629491 2.138024 32.550000 33.930000 34.250000 37.710000 41.950000
housing_median_age 16340.0 28.663525 12.556964 1.000000 18.000000 29.000000 37.000000 52.000000
total_rooms 16340.0 2622.775949 2169.626947 11.000000 1449.000000 2119.500000 3121.250000 39320.000000
total_bedrooms 16340.0 535.390208 418.478468 3.000000 296.000000 433.000000 644.000000 6445.000000
population 16340.0 1415.480171 1101.940296 3.000000 785.000000 1164.000000 1711.000000 28566.000000
households 16340.0 496.932313 379.070022 3.000000 280.000000 408.000000 600.250000 6082.000000
median_income 16340.0 3.869024 1.901773 0.499900 2.555600 3.533000 4.744075 15.000100
median_house_value 16340.0 206921.906977 115662.820190 14999.000000 119300.000000 179700.000000 265900.000000 500001.000000
people_per_household 16340.0 2.946757 1.291597 0.692308 2.428869 2.816514 3.282433 83.171429
rooms_per_household 16340.0 5.413407 2.023690 0.846154 4.443662 5.229614 6.050922 52.848214
bedrooms_per_people 16340.0 0.402226 0.202989 0.011222 0.315537 0.372646 0.443428 7.925926
bedrooms_per_rooms 16340.0 0.213188 0.057896 0.100000 0.175588 0.203260 0.239666 1.000000

Analyse kategorialer Variablen

# summary statistics for all categorical columns
train_dataset.describe(include=['category']).transpose()
count unique top freq
ocean_proximity 16340 5 <1H OCEAN 7216

Die Variable ocean_proximity besitzt 5 Ausprägungen, wovon die Ausprägung <1H OCEAN mit einer Anzhl von 7215 am häufigsten vorkommt.

train_dataset['ocean_proximity'].value_counts()
<1H OCEAN     7216
INLAND        5216
NEAR OCEAN    2107
NEAR BAY      1796
ISLAND           5
Name: ocean_proximity, dtype: int64

Hier wird die Häufigkeit einzelner Ausprägungen von ocean_proximity dargestellt. Die Ausprägung ISLAND kommt über das gesamte Datenset mit über 20000 Datensätzen nur fünfmal vor. In Relation zur Gesamtanzahl der Datensätzen ist dies sehr wenig, wodurch die Aussagekraft dieser Ausprägung vermutlich eher gering ist. Daher werden die Datensätze mit Ausprägung ISLAND mit der Ausprägung NEAR OCEAN vereint.

train_dataset.loc[train_dataset.ocean_proximity == 'ISLAND']
longitude latitude housing_median_age total_rooms total_bedrooms population households median_income median_house_value ocean_proximity people_per_household rooms_per_household bedrooms_per_people bedrooms_per_rooms
8318 -118.48 33.43 29 716 214 422 173 2.6042 287500.0 ISLAND 2.439306 4.138728 0.507109 0.298883
8316 -118.32 33.33 52 2127 512 733 288 3.3906 300000.0 ISLAND 2.545139 7.385417 0.698499 0.240715
8315 -118.33 33.34 52 2359 591 1100 431 2.8333 414700.0 ISLAND 2.552204 5.473318 0.537273 0.250530
8314 -118.32 33.35 27 1675 521 744 331 2.1579 450000.0 ISLAND 2.247734 5.060423 0.700269 0.311045
8317 -118.32 33.34 52 996 264 341 160 2.7361 450000.0 ISLAND 2.131250 6.225000 0.774194 0.265060
# write out a dict with the mapping of old to new
remap_cat_dict = {
    'NEAR OCEAN': 'NEAR OCEAN',
    'NEAR BAY': 'NEAR BAY',
    'ISLAND': 'NEAR OCEAN',
    'INLAND': 'INLAND',
    '<1H OCEAN': '<1H OCEAN' }

train_dataset.ocean_proximity = train_dataset.ocean_proximity.map(remap_cat_dict).astype('category')

Exploratory Data Anaysis

Dieses Kapitel dient der ausführlichen Datenanalyse der Testdaten. Hierfür werden vor allem geeignete Visualisierungen in Form von Diagrammen verwendet.

erster Überblick

# summary statistics for all numerical columns
round(train_dataset.describe(),2).transpose()
count mean std min 25% 50% 75% max
longitude 16340.0 -119.56 2.00 -124.35 -121.79 -118.49 -118.00 -114.47
latitude 16340.0 35.63 2.14 32.55 33.93 34.25 37.71 41.95
housing_median_age 16340.0 28.66 12.56 1.00 18.00 29.00 37.00 52.00
total_rooms 16340.0 2622.78 2169.63 11.00 1449.00 2119.50 3121.25 39320.00
total_bedrooms 16340.0 535.39 418.48 3.00 296.00 433.00 644.00 6445.00
population 16340.0 1415.48 1101.94 3.00 785.00 1164.00 1711.00 28566.00
households 16340.0 496.93 379.07 3.00 280.00 408.00 600.25 6082.00
median_income 16340.0 3.87 1.90 0.50 2.56 3.53 4.74 15.00
median_house_value 16340.0 206921.91 115662.82 14999.00 119300.00 179700.00 265900.00 500001.00
people_per_household 16340.0 2.95 1.29 0.69 2.43 2.82 3.28 83.17
rooms_per_household 16340.0 5.41 2.02 0.85 4.44 5.23 6.05 52.85
bedrooms_per_people 16340.0 0.40 0.20 0.01 0.32 0.37 0.44 7.93
bedrooms_per_rooms 16340.0 0.21 0.06 0.10 0.18 0.20 0.24 1.00
import plotly.graph_objects as go

fig = go.Figure(data=go.Scattergeo(
        lon = train_dataset['longitude'],
        lat = train_dataset['latitude'],
        text = 'median house value in district: ' + train_dataset['median_house_value'].astype(str) + '$',
        marker_color = train_dataset['median_house_value'], marker_size=2, marker_colorscale = 'Inferno', marker_colorbar= {'title': 'Median House Value'}
        ))

fig.update_layout(
        title = 'Median House Value per District in California',
        geo_scope='usa',
    )
fig.show()

Die Datensätze stammen alle aus Kalifornien, USA.
Die Werte Longitude und Latitude unterscheiden sich daher kaum, wie sich auch bereits im Kapitel deskriptive Statistik gezeigt hat. Daher sind sie als Feature eher ungeeignet.
Allerdings scheint die Entfernung zum Ozean eine wichtige Rolle im Bezug auf den median_house_value zu spielen. Ocean_proximity scheint daher als feature interessant und sollte genauer untersucht werden.

train_dataset = train_dataset.drop(['longitude', 'latitude'], axis = 1)
sns.pairplot(train_dataset[['housing_median_age', 'total_rooms', 'total_bedrooms', 'population', 'households', 'median_house_value']])
<seaborn.axisgrid.PairGrid at 0x1ec11b17490>
Error in callback <function flush_figures at 0x000001EC10073040> (for post_execute):
---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
C:\ProgramData\Anaconda3\lib\site-packages\matplotlib_inline\backend_inline.py in flush_figures()
    119         # ignore the tracking, just draw and close all figures
    120         try:
--> 121             return show(True)
    122         except Exception as e:
    123             # safely show traceback if in IPython, else raise

C:\ProgramData\Anaconda3\lib\site-packages\matplotlib_inline\backend_inline.py in show(close, block)
     39     try:
     40         for figure_manager in Gcf.get_all_fig_managers():
---> 41             display(
     42                 figure_manager.canvas.figure,
     43                 metadata=_fetch_figure_metadata(figure_manager.canvas.figure)

C:\ProgramData\Anaconda3\lib\site-packages\IPython\core\display.py in display(include, exclude, metadata, transient, display_id, *objs, **kwargs)
    318             publish_display_data(data=obj, metadata=metadata, **kwargs)
    319         else:
--> 320             format_dict, md_dict = format(obj, include=include, exclude=exclude)
    321             if not format_dict:
    322                 # nothing to display (e.g. _ipython_display_ took over)

C:\ProgramData\Anaconda3\lib\site-packages\IPython\core\formatters.py in format(self, obj, include, exclude)
    178             md = None
    179             try:
--> 180                 data = formatter(obj)
    181             except:
    182                 # FIXME: log the exception

C:\ProgramData\Anaconda3\lib\site-packages\decorator.py in fun(*args, **kw)
    230             if not kwsyntax:
    231                 args, kw = fix(args, kw, sig)
--> 232             return caller(func, *(extras + args), **kw)
    233     fun.__name__ = func.__name__
    234     fun.__doc__ = func.__doc__

C:\ProgramData\Anaconda3\lib\site-packages\IPython\core\formatters.py in catch_format_error(method, self, *args, **kwargs)
    222     """show traceback on failed format call"""
    223     try:
--> 224         r = method(self, *args, **kwargs)
    225     except NotImplementedError:
    226         # don't warn on NotImplementedErrors

C:\ProgramData\Anaconda3\lib\site-packages\IPython\core\formatters.py in __call__(self, obj)
    339                 pass
    340             else:
--> 341                 return printer(obj)
    342             # Finally look for special method names
    343             method = get_real_method(obj, self.print_method)

C:\ProgramData\Anaconda3\lib\site-packages\IPython\core\pylabtools.py in print_figure(fig, fmt, bbox_inches, base64, **kwargs)
    149         FigureCanvasBase(fig)
    150 
--> 151     fig.canvas.print_figure(bytes_io, **kw)
    152     data = bytes_io.getvalue()
    153     if fmt == 'svg':

C:\ProgramData\Anaconda3\lib\site-packages\matplotlib\backend_bases.py in print_figure(self, filename, dpi, facecolor, edgecolor, orientation, format, bbox_inches, pad_inches, bbox_extra_artists, backend, **kwargs)
   2292             if bbox_inches:
   2293                 if bbox_inches == "tight":
-> 2294                     bbox_inches = self.figure.get_tightbbox(
   2295                         renderer, bbox_extra_artists=bbox_extra_artists)
   2296                     if pad_inches is None:

C:\ProgramData\Anaconda3\lib\site-packages\matplotlib\figure.py in get_tightbbox(self, renderer, bbox_extra_artists)
   1626 
   1627         for a in artists:
-> 1628             bbox = a.get_tightbbox(renderer)
   1629             if bbox is not None and (bbox.width != 0 or bbox.height != 0):
   1630                 bb.append(bbox)

C:\ProgramData\Anaconda3\lib\site-packages\matplotlib\artist.py in get_tightbbox(self, renderer)
    353             The enclosing bounding box (in figure pixel coordinates).
    354         """
--> 355         bbox = self.get_window_extent(renderer)
    356         if self.get_clip_on():
    357             clip_box = self.get_clip_box()

C:\ProgramData\Anaconda3\lib\site-packages\matplotlib\patches.py in get_window_extent(self, renderer)
    621 
    622     def get_window_extent(self, renderer=None):
--> 623         return self.get_path().get_extents(self.get_transform())
    624 
    625     def _convert_xy_units(self, xy):

C:\ProgramData\Anaconda3\lib\site-packages\matplotlib\path.py in get_extents(self, transform, **kwargs)
    612         from .transforms import Bbox
    613         if transform is not None:
--> 614             self = transform.transform_path(self)
    615         if self.codes is None:
    616             xys = self.vertices

C:\ProgramData\Anaconda3\lib\site-packages\matplotlib\transforms.py in transform_path(self, path)
   1606         that began as line segments.
   1607         """
-> 1608         return self.transform_path_affine(self.transform_path_non_affine(path))
   1609 
   1610     def transform_path_affine(self, path):

C:\ProgramData\Anaconda3\lib\site-packages\matplotlib\transforms.py in transform_path_affine(self, path)
   1616         ``transform_path_affine(transform_path_non_affine(values))``.
   1617         """
-> 1618         return self.get_affine().transform_path_affine(path)
   1619 
   1620     def transform_path_non_affine(self, path):

C:\ProgramData\Anaconda3\lib\site-packages\matplotlib\transforms.py in get_affine(self)
   2444             return self._b.get_affine()
   2445         else:
-> 2446             return Affine2D(np.dot(self._b.get_affine().get_matrix(),
   2447                                    self._a.get_affine().get_matrix()))
   2448 

C:\ProgramData\Anaconda3\lib\site-packages\matplotlib\transforms.py in __init__(self, matrix, **kwargs)
   1912             # A bit faster than np.identity(3).
   1913             matrix = IdentityTransform._mtx.copy()
-> 1914         self._mtx = matrix.copy()
   1915         self._invalid = 0
   1916 

KeyboardInterrupt: 

Die hier abgebildeten, unabhängigen Variablen weisen eine starke Korrelation untereinander auf (bis auf housing_median_age).
In Bezug auf die abhängige Variable median_house_value zeigen die Plots zunächst keinen starken Zusammenhang.
Interessant sieht das Histogramme von median_house_value aus, da es am äußeren rechten Rand nochmal hohe Häufigkeiten haben. Dieses Histogramme wird später nochmal genauer analysisiert.

sns.pairplot(train_dataset[['median_income', 'people_per_household', 'rooms_per_household', 'bedrooms_per_people', 'bedrooms_per_rooms', 'median_house_value']]);
_images/regression_57_0.png

Zwischen median_income und median_house_value ist ein Zusammenhang erkennbar. Auch das Histogramm von median_income zeigt eine hohe Säule am rechten Rand.
Auch für die selbst definierten Variablen ist es schwierig, ein linearen Zusammenhang zu median_house_value zu erkennen.

fig = sns.displot(data= train_dataset, x='median_house_value')
_images/regression_59_0.png

Die Verteilung von Median_House_value nicht normalverteilt, rechts im hohen Preissegment ist die Häufigkeit nochmal extrem hoch. Über 500000\( exisiteren keine weiteren Hauswerte. Möglicherweise wurde der Datensatz mit Werten über 500000\) abgeschnitten und für alle höheren Werte, der Wert 500000$ gesetzt.

fig = sns.displot(data= train_dataset, x='median_income')
_images/regression_61_0.png

Womöglich wurde auch die Variable median_income nach oben abgeschnitten.

Korrelationsanalyse

# Create correlation matrix for numerical variables
corr_matrix = train_dataset.corr()
corr_matrix['median_house_value'].sort_values(ascending=False)
median_house_value      1.000000
median_income           0.692088
rooms_per_household     0.176506
total_rooms             0.130850
housing_median_age      0.105574
bedrooms_per_people     0.077531
households              0.063122
total_bedrooms          0.047633
population             -0.026782
longitude              -0.043011
latitude               -0.147308
people_per_household   -0.157646
bedrooms_per_rooms     -0.253006
Name: median_house_value, dtype: float64
# Make a pretty heatmap

# Use a mask to plot only part of a matrix
mask = np.zeros_like(corr_matrix)
mask[np.triu_indices_from(mask)]= True

# Change size
plt.subplots(figsize=(11, 15))

# Build heatmap with additional options
heatmap = sns.heatmap(corr_matrix, 
                      mask = mask, 
                      square = True, 
                      linewidths = .5,
                      cmap = 'coolwarm',
                      cbar_kws = {'shrink': .6,
                                'ticks' : [-1, -.5, 0, 0.5, 1]},
                      vmin = -1,
                      vmax = 1,
                      annot = True,
                      annot_kws = {"size": 10})
_images/regression_65_0.png

Laut der Korrelationsmatrix einigen sich vor allem folgende Variablen:
median_income weist mit einem Wert von 0,69 die höchste Korrelation auf, danach folgen bedrooms_per_rooms mit einer negativen Korrelation in Höhe von -0,25 und anschließend people_per_household ebenfalls mit einer negativen Korrelation in Höhe von -0,16.
Danach folgt rooms_per_household, allerdings kann diese Variable nicht in Kombination mit den beiden zuvor genannten genutzt werden, da sie sich aus total_rooms und households zusammensetzt. Diese Werte werden jedoch bereits in bedrooms_per_rooms und people_per_household verwendet, daher ist die zusätzliche Verwendung von rooms_per_household nicht sinnvoll.
housing_median_age könnte mit einer Korrelation von 0,11 noch einen positiven Einfluss auf das Modell haben.

Analyse der kategorialen Variable Ocean_Proximity

sns.boxenplot(data=train_dataset, x='ocean_proximity', y='median_house_value')
<AxesSubplot:xlabel='ocean_proximity', ylabel='median_house_value'>
_images/regression_68_1.png
sns.stripplot(data=train_dataset, x='ocean_proximity', y='median_house_value', size=2)
<AxesSubplot:xlabel='ocean_proximity', ylabel='median_house_value'>
_images/regression_69_1.png

Die Ausprägungen NEAR BAY und NEAR OCEAN unterscheiden sich hinsichtlich median_house_value kaum. Daher können diese beiden Ausprägungen zusammengefasst werden.

sns.displot(data=train_dataset, x='median_house_value', hue='ocean_proximity', kind='kde')
<seaborn.axisgrid.FacetGrid at 0x230d8270f40>
_images/regression_71_1.png

Auch in der Häufigkeitsverteilung zeigen NEAR BAY und NEAR OCEAN kaum Unterschiede.

# write out a dict with the mapping of old to new
remap_cat_dict = {
    'NEAR OCEAN': 'NEAR COAST',
    'NEAR BAY': 'NEAR COAST',
    'INLAND': 'INLAND',
    '<1H OCEAN': '<1H OCEAN' }

train_dataset.ocean_proximity = train_dataset.ocean_proximity.map(remap_cat_dict).astype('category')
sns.boxplot(data=train_dataset, x='ocean_proximity', y='median_house_value')
<AxesSubplot:xlabel='ocean_proximity', ylabel='median_house_value'>
_images/regression_74_1.png

Auch NEAR COAST und <1H OCEAN zeigen eine ähnliche Verteilung des Medians und der Quartile. Das 3. Quartil von NEAR COAST zieht sich jedoch noch um ca. 50000$ weiter nach oben, als das von <1H OCEAN, daher bleibt die Differenzierung dieser beiden Ausprägungen noch bestehen.
Bei der Ausprägung INLAND sind viele Ausreiser nach oben zu erkennen.

Analyse numerischer Variablen

# check relationship with a joint plot
sns.jointplot(x="median_income", y="median_house_value", hue="ocean_proximity", data=train_dataset);
_images/regression_77_0.png
# check relationship with a joint plot
sns.lmplot(x="median_income", y="median_house_value", data=train_dataset);
_images/regression_78_0.png

Wie bereits aus der Korrelationsmatrix zu schließen, ist ein linearer Zusammenhang zwischen median_income und median_house_value erkennbar. Die “Linie” mit Datenpunkten bei 500000\( ergibt sich aus der anzunehmenden Begrenzung der median_house_value_ auf 500000\).
Auch innerhalb der einzelenden Ausprägungen von ocean_proximity ist eine lineare Korrelation von median_income und median_house_value erkennbar.

sns.relplot(data=train_dataset, x='total_rooms', y='median_house_value', hue='ocean_proximity', col='ocean_proximity')
<seaborn.axisgrid.FacetGrid at 0x193868f47c0>
_images/regression_80_1.png

Ein linearer Zusammenhang von total_rooms zu median_house_value ist nicht wirklich zu erkennen. Die Variable ist daher für eine lineare Regression eher ungeeignet.

# Plot with Plotly Express
px.scatter(x=train_dataset['bedrooms_per_rooms'], y=train_dataset['median_house_value'], opacity=0.65, 
                trendline='ols', trendline_color_override='darkred')
px.scatter(x=train_dataset['people_per_household'], y=train_dataset['median_house_value'], opacity=0.65, 
                trendline='ols', trendline_color_override='darkred')
sns.lmplot(data=train_dataset, x='people_per_household', y='median_house_value')
<seaborn.axisgrid.FacetGrid at 0x1a1836ebca0>
_images/regression_84_1.png

bedrooms_per_rooms und people_per_household zeigen einen negativen linearen Zusammenhang. Dies ists gut an der Linie im Plot verdeutlicht erkennbar. Bereits in der Korrelationsmatrix war der negative Zusammenhang dieser beiden Variablen zu median_house_value sichtbar. Die Variablen werden nachfolgend in die linearen Regressionsmodelle aufgenommen.

Modellierung

Statsmodel

Model 1 - OLS Regression

Als erstes Modell wird statsmodels OLS Regression verwendet.
Die sogenannten Filter Methode ist eine von drei Möglichkeiten bei der Feature Selektion. Gemäß dieser werden beispielsweise auf Basis einer Korrelationsmatrix die besten Feature ausgewählt und diese dem Modell zur Verfügung gestelltn.
Folgende Features wurden als geeignet identifiziert: median_income, ocean_proximity, bedrooms_per_rooms, people_per_household und eventuell housing_median_age.
Im ersten Modell werden die vier zuerst genannten Features verwendet.

# Fit Model
lm = smf.ols(formula='median_house_value ~ median_income + ocean_proximity + bedrooms_per_rooms + people_per_household' , data=train_dataset).fit()
# Full summary
lm.summary()
OLS Regression Results
Dep. Variable: median_house_value R-squared: 0.619
Model: OLS Adj. R-squared: 0.619
Method: Least Squares F-statistic: 5312.
Date: Tue, 18 Jan 2022 Prob (F-statistic): 0.00
Time: 01:03:32 Log-Likelihood: -2.0580e+05
No. Observations: 16340 AIC: 4.116e+05
Df Residuals: 16334 BIC: 4.116e+05
Df Model: 5
Covariance Type: nonrobust
coef std err t P>|t| [0.025 0.975]
Intercept 2.285e+04 4519.014 5.056 0.000 1.4e+04 3.17e+04
ocean_proximity[T.INLAND] -6.729e+04 1415.385 -47.545 0.000 -7.01e+04 -6.45e+04
ocean_proximity[T.NEAR COAST] 1.676e+04 1428.672 11.729 0.000 1.4e+04 1.96e+04
median_income 4.338e+04 406.094 106.829 0.000 4.26e+04 4.42e+04
bedrooms_per_rooms 3.031e+05 1.3e+04 23.248 0.000 2.78e+05 3.29e+05
people_per_household -1.049e+04 435.520 -24.097 0.000 -1.13e+04 -9640.930
Omnibus: 4763.697 Durbin-Watson: 1.971
Prob(Omnibus): 0.000 Jarque-Bera (JB): 22109.483
Skew: 1.349 Prob(JB): 0.00
Kurtosis: 8.019 Cond. No. 129.


Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
lm.params
Intercept                         22847.458850
ocean_proximity[T.INLAND]        -67294.434300
ocean_proximity[T.NEAR COAST]     16757.013965
median_income                     43382.483758
bedrooms_per_rooms               303138.303742
people_per_household             -10494.596699
dtype: float64

Die summary() Methode berechnet die wichtigsten Kennzahlen des Modells.
R-squared ist eine Kennzahl zur Beurteilung der Anpassungsgüte eines Modells und nimmt einen Wert zwischen 0% und 100% an. Als Bezugsbasis wird der Durchschnitt verwendet.
lm erreicht einen R-squared Wert in Höhe von 63,2%. Somit können 63,2% der median_house_value Streuung
In statistischen Methoden wird meist lieber der adjusted R-squared genutzt. Da dieser die degrees of freedom miteinbezieht, lässt sich durch adj. R-squared einen besseren Rückschluss auf die Gesamtpopulation der Datensätze ziehen. Dieser ist meistens etwas schlechter als R-squared. In diesem Modell sind beide Werte gleich hoch.
Mithilfe der Kennzahlen AIC und BIC soll ein overfitting des Modells verhindert werden. Für jedes zusätzliche Feature wird eine “Bestrafung” in die Kennzahlen miteinberechnet. Je kleiner die Werte, desto besser. In diesem Modell sind die Werte jedoch relativ hoch.
F-statistics ist die Menge der systematischen Varianz (MSM) geteilt durch die Menge der unsystematischen Varianz (MSR). Die Kennzahl gibt an, in welcher Höhe das Modell die Ergebnisausgabe der abhängigen Variable verbessert hat, verglichen zu Ungenauigkeit im Modell. Je höher F-statistics, desto besser. Ein F-statistics Wert in Höhe von 5618, wie in Modell lm, ist schon sehr gut.
Der Y-Achsenschnittpunkt des Modells liegt bei 61394,65. Die Steigung der einzelnen Features ist in der Spalte coef abzulesen. Zu interpretieren sind diese folgendermaßen: Vergrößert sich der Wert in Variable median_income um 1 und alle anderen Werte bleiben konstant, so ist das Ergebnis um 42189,93 größer als zuvor.
Die p-values der Features liegen alle bei 0,00. Die Nullhypothese (ß=0 => Feature hat keinen Einfluss auf die abhängige Variable) kann demnach für alle Features verworfen werden.
Die p-values für Jarque-Bera und Omnibus liegen ebenfalls bei 0,00. Die Nullhypothese, die Fehler seien normalverteilt, wird in diesem Falle verworfen. Es wird also angenommen, die Residuals seien nicht normalverteilt. Bei einer Datensetgröße > 50, kann dies aufgrund des Central Limit Theorems jedoch vernachlässigt werden.
Der Durbin-Watson Wert liegt bei knapp 2. Somit kann angenommen werden, dass die Fehler des Modells keine Korrelation aufweisen.

Regression Diagnostics
fig = sm.graphics.influence_plot(lm, criterion="cooks")
fig.tight_layout(pad=1.0)
_images/regression_94_0.png

Mithilfe des Influence Plots werden mögliche Ausreiser als blaue Punkte in der Größe entsprechend ihres Einflusses identifiziert. Dies können Werte mit hoher Abweichung in ihrem Y-Wert (outliers) oder X-Wert (high leverage) sein. In einer linearen Regression können solche Ausreise einen großen Einfluss auf das Gesamtmodell haben und somit die Performance des gesamten Modells erheblich verschlechtern. Als Faustregel werden Datenpunkte mit einer Cooks Distance in Höhe von 4/n (wobei n die Anzahl der Datensätze ist) als mögliche Ausreiser identifiziert.

# obtain Cook's distance 
lm_cooksd = lm.get_influence().cooks_distance[0]

# get length of df to obtain n
n = len(train_dataset["median_income"])
print('n = ', n)

# calculate critical d
critical_d = 4/n
print('Critical Cooks distance:', critical_d)

# identification of potential outliers with leverage
out_d = lm_cooksd > critical_d

# output potential outliers with leverage
print(train_dataset.index[out_d], "\n", 
    lm_cooksd[out_d])
n =  16340
Critical Cooks distance: 0.0002447980416156671
Int64Index([15778,  6776, 11511,  7263, 15723, 14756, 15687,  4638, 14806,
            15815,
            ...
            15809,  8866,   199, 18498,    59, 10749,  8857, 13154,  5470,
            16642],
           dtype='int64', length=921) 
 [1.48759293e-03 2.66836227e-04 4.78884475e-04 5.45180812e-04
 5.82637029e-04 5.55166261e-03 5.93629188e-04 9.64446395e-04
 7.14112179e-04 8.23043391e-04 1.37179384e-03 9.31198962e-04
 5.28353433e-04 3.14064566e-04 3.70703445e-04 7.92176597e-04
 5.76686718e-04 2.53437155e-04 5.35087643e-04 4.34543876e-04
 2.56743264e-04 4.37260421e-04 3.15946778e-04 3.85848620e-04
 2.78654618e-04 9.85584673e-04 4.69834930e-02 5.19188525e-04
 4.33272705e-04 2.46198971e-04 3.62984273e-04 3.72530337e-04
 2.66246300e-04 8.19396201e-04 4.87217083e-04 1.63877589e-03
 2.18039096e-03 1.10631997e-03 5.10220751e-04 2.98320853e-04
 4.98424061e-04 2.83831187e-03 4.23210697e-04 4.84462529e-04
 8.35832817e-04 2.56077726e-03 6.60015026e-04 3.00979233e-04
 3.89758021e-04 9.03255145e-04 3.76539069e-04 4.30606581e-04
 3.28854120e-04 3.86905704e-04 7.61360057e-04 3.48755331e-04
 4.90570754e-04 5.32416641e-04 4.53912402e-04 5.91744313e-04
 3.61448425e-04 4.27661004e-04 4.01913126e-04 7.92934770e-04
 4.08537643e-04 4.60726908e-04 1.17463913e-03 2.51060698e-04
 4.80930856e-04 5.12470023e-04 5.65809255e-04 3.27401509e-04
 6.85242374e-04 1.21787762e-03 7.24325239e-04 2.92813642e-04
 3.73003816e-04 4.68626337e-04 1.16110752e-03 4.63062325e-04
 4.25016890e-04 4.83688775e-04 5.35658239e-04 1.94858180e-03
 2.89813290e-04 8.85370793e-04 7.17970203e-04 6.42891705e-04
 3.14955277e-04 1.02215153e-03 3.56207358e-04 8.67510310e-04
 1.10491710e-03 7.99219333e-04 2.62696017e-04 2.88666981e-04
 5.07854509e-04 8.83931891e-04 2.45081661e-04 1.06449798e-03
 3.70938773e-04 3.09018081e-04 6.57379410e-04 1.06179504e-03
 3.02894187e-04 9.26555582e-04 1.01293834e-03 3.20886325e-04
 2.31556876e-03 3.21853302e-04 3.36049381e-04 2.50133964e-04
 5.13847455e-04 2.52642306e-04 3.13725249e-04 7.12813234e-04
 4.91857937e-04 5.77017000e-04 6.85531264e-04 2.60496718e-04
 2.85522112e-04 2.61791680e-04 3.35125004e-04 7.60609675e-04
 5.47063429e-04 8.90184001e-04 9.82453636e-04 3.85061491e-04
 5.98611688e-04 4.21507557e-04 5.80036377e-04 2.97849695e-04
 6.48910258e-04 4.04532613e-04 3.12791855e-02 7.55680452e-04
 3.20510714e-04 2.91560940e-04 3.88296383e-04 2.61226212e-04
 7.36742340e-04 7.57230958e-04 6.60575314e-04 5.89682358e-04
 2.53674156e-04 5.67379890e-04 4.75221026e-04 3.19645157e-04
 2.78539449e-04 3.63618083e-04 2.74057139e-04 5.54655710e-04
 2.54335893e-04 2.75933612e-04 4.70547153e-04 6.15449463e-04
 2.62179019e-04 3.00882646e-04 4.13638589e-04 2.65252379e-04
 9.58798052e-04 3.59247990e-04 1.25803553e-03 4.35828212e-03
 2.77124480e-04 3.30386011e-04 6.10246492e-04 3.95647548e-04
 3.18201721e-04 3.56188388e-04 2.78201241e-04 2.45863813e-04
 5.02007922e-04 2.50010455e-04 6.14141780e-04 2.92218505e-04
 2.53697899e-04 2.46267670e-04 2.67098403e-04 7.51479956e-04
 7.81424778e-04 3.27290114e-04 8.80851341e-04 1.20674591e-03
 3.03575822e-03 5.18664791e-04 3.89595116e-04 2.68812640e-04
 2.59722638e-04 1.20517420e-03 6.24208641e-04 3.72779920e-03
 8.51415384e-04 3.72957326e-04 3.05334929e-04 5.52136850e-04
 2.98068484e-03 3.29411537e-04 3.40411907e-04 2.56222465e-04
 6.31335952e-04 3.11886294e-03 2.61795401e-04 3.49696722e-04
 2.74889698e-04 4.24423862e-04 5.33480402e-03 3.76297264e-03
 2.83074019e-03 8.29506869e-04 1.36747265e-03 4.79332644e-04
 4.70922385e-04 3.90781122e-04 2.64913644e-04 5.65879731e-04
 5.98086264e-04 2.94669543e-04 2.63610289e-04 3.68587927e-04
 3.99298783e-04 4.35823599e-04 9.11267981e-04 9.39988920e+00
 3.39861014e-04 2.69052780e-04 3.61433050e-04 6.04443335e-04
 2.66277640e-04 7.66984878e-04 3.10471732e-04 4.34022341e-02
 5.09218138e-04 3.68185782e-04 2.18766380e-03 3.56981125e-04
 4.31902936e-04 1.39521400e-03 3.20738655e-04 4.98230233e-04
 4.77448218e-04 3.32692672e-04 2.94224674e-02 9.53393326e-04
 4.06615406e-04 5.79773700e-04 3.51514909e-04 3.61553795e-04
 1.22957130e-03 4.63021859e-04 2.69086838e-04 4.12697044e-04
 3.51813451e-04 9.70081562e-04 2.55405842e-04 1.29836196e-03
 3.14407763e-03 4.27009366e-04 8.02132198e-04 2.68103122e-04
 2.88662584e-04 5.24588063e-04 2.47282051e-04 3.09316188e-04
 9.85607781e-04 3.62042318e-04 2.69988391e-04 2.60804951e-04
 2.95435073e-04 2.64265161e-03 2.73619549e-04 3.09846370e-04
 1.08953714e-03 3.19858062e-03 3.61945601e-04 4.48200055e-04
 4.46775381e-04 2.45696522e-04 3.29958181e-04 4.37556828e-04
 5.94720981e-04 8.15616897e-04 2.94586664e-04 2.46671711e-04
 3.67944042e-04 5.20079922e-04 2.61685051e-04 3.00405843e-03
 2.95447880e-04 2.69710091e-04 5.84435890e-04 7.01102896e-04
 3.06679338e-03 5.43870147e-03 3.74467022e-04 2.94450188e-04
 2.53604644e-04 3.86371760e-04 3.37832122e-04 4.08310873e-04
 2.57857913e-04 2.53449617e-04 4.35238135e-04 9.32109329e-04
 4.19634636e-04 3.70736472e-03 6.60444604e-04 2.13403879e-03
 4.47318394e-04 3.93509961e-03 1.71125362e-03 3.14016243e-04
 8.64171619e-04 3.54606811e-04 3.03053219e-04 4.38959252e-04
 1.86291326e-03 4.18578111e-04 6.23666239e-04 2.74507153e-04
 5.55770133e-04 2.39766662e-03 5.29707607e-04 1.47048722e-03
 5.72050204e-04 5.91643897e-03 6.45921796e-04 1.71379171e-03
 3.29776306e-04 4.08631064e-04 1.03297005e-03 5.77032046e-04
 2.76044202e-04 6.17879439e-04 4.20098365e-04 1.41320439e-03
 2.63198371e-04 2.39553721e-03 3.56444674e-04 3.76120928e-04
 3.75049105e-04 4.20190411e-03 3.52751919e-04 8.64347341e-04
 2.14599731e-03 5.21634440e-04 5.83742875e-04 2.69107084e-04
 5.94255954e-04 6.08876932e-04 8.17456580e-04 4.13993717e-04
 3.67213098e-04 7.19449438e-04 2.52664572e-04 3.51423970e-04
 3.71139725e-04 4.30278077e-04 1.10036644e-03 4.48525821e-04
 2.97875784e-04 2.53128372e-04 6.71543161e-04 2.57889557e-04
 4.37631343e-04 5.31753783e-04 4.36362942e-04 3.48171741e-04
 5.33669610e-03 6.27844689e-04 3.30353881e-04 2.63316932e-04
 1.93430902e-03 3.84002969e-04 3.10408934e-04 8.02639776e-04
 3.60331666e-04 3.55062676e-04 2.93513559e-04 7.77477806e-04
 2.78901271e-04 2.56911171e-04 2.62926487e-04 2.66958932e-04
 3.18433088e-04 4.64211041e-04 3.30739046e-04 2.52024154e-04
 1.27278215e-03 2.84426835e-04 7.04396718e-04 6.33891174e-04
 3.39793952e-04 2.99236648e-04 4.45190246e-04 5.35835774e-04
 4.22579904e-04 1.14102901e-03 1.00360134e-03 3.28224562e-04
 2.79656610e-04 3.27410721e-04 7.39804588e-04 3.89886247e-04
 3.33262252e-04 3.94098964e-04 4.18819969e-04 8.19442679e-04
 3.89986966e-04 5.36623487e-04 3.96495856e-04 2.53190918e-03
 7.24780554e-04 2.76551422e-03 2.68771101e-04 3.54448329e-04
 1.13845602e-03 2.46529930e-04 4.32895348e-03 4.15819743e-04
 3.93225300e-04 3.82160961e-03 5.93387422e-04 6.77249628e-04
 4.68111655e-04 1.48386533e-03 2.70829420e-04 2.60996883e-03
 3.80995269e-04 5.33430607e-04 5.13911021e-04 2.56822250e-04
 9.00464300e-04 5.89936555e-04 1.00222242e-03 4.49752787e-04
 4.85344597e-03 4.42173403e-04 3.74656557e-04 4.05218668e-04
 2.49437119e-04 8.48112319e-04 9.52066468e-04 3.17766109e-04
 3.21040279e-04 2.67941768e-04 2.16147046e-03 5.23891516e-04
 3.89009713e-03 5.82745086e-04 6.78943360e-04 3.07943494e-04
 5.01012970e-04 4.48561295e-04 1.28186822e-03 8.61789035e-04
 4.34727316e-04 2.73804253e-04 1.42202399e-01 3.50860752e-04
 2.77507369e-03 4.14111841e-04 4.34040630e-04 8.47409581e-04
 4.73707382e-04 9.36056215e-04 4.39638949e-04 2.47091843e-04
 2.94815965e-03 3.41359419e-04 2.58914669e-04 1.29921590e-03
 3.43443742e-04 7.08910616e-04 2.51070091e-04 7.51519537e-04
 3.61757210e-04 2.28697927e-03 3.46303378e-04 2.72064007e-04
 4.02842045e-04 1.08821274e-03 2.48982171e-03 2.53807767e-04
 2.76891853e-04 2.47227261e-04 3.45705294e-04 4.99565882e-04
 2.35401915e-03 2.70401709e-04 1.09188172e-02 2.91184063e-04
 5.50205439e-04 3.35361471e-04 7.99127506e-04 4.69299374e-04
 2.63469833e-04 7.40940278e-04 2.64031971e-04 3.29470784e-04
 1.47023459e-03 3.40530222e-04 4.65887402e-03 4.50849207e-04
 3.52524725e-04 8.97162356e-03 5.70524019e-04 7.51568069e-04
 2.50790326e-04 2.85441357e-03 6.78405443e-03 3.43573733e-04
 3.36485914e-04 2.98686005e-04 8.58965733e-04 2.54924272e-04
 8.30030938e-04 2.56453063e-03 3.93964524e-04 3.95499423e-04
 3.79750484e-04 5.50201808e-04 1.71339260e-03 3.37091835e-04
 4.40192757e-04 8.24478158e-04 3.12070719e-04 1.35764062e-03
 3.43487813e-04 3.12996792e-04 5.68496869e-04 3.98348038e-04
 4.46630274e-04 4.95629626e-04 2.57410000e-04 8.08578807e-04
 2.64273317e-04 2.82160955e-04 3.65064309e-04 2.70263308e-04
 2.45828059e-04 8.38319883e-04 3.25170223e-04 3.06324121e-04
 3.33351819e-04 6.77037097e-03 2.75609130e-04 1.01649062e-03
 5.04145849e-04 5.15573187e-04 3.88917919e-04 1.46058927e-03
 3.71213447e-04 1.08778812e-03 2.56939691e-04 6.19395605e-04
 5.94834315e-04 4.06971618e-04 2.78217893e-03 3.23101550e-04
 7.74127917e-04 1.76364447e-03 8.36296107e-04 3.80469230e-04
 2.70007674e-04 1.17371929e-03 6.95576355e-04 3.32710828e-04
 2.59662126e-04 5.69312329e-04 9.26900479e-04 8.42265531e-04
 3.10165863e-04 2.70328570e-03 7.71394059e-04 2.68445474e-04
 3.11973612e-04 3.32461369e-04 9.80372135e-04 8.81917284e-04
 7.55373959e-04 3.29644215e-04 6.21269904e-04 3.76872079e-04
 3.50203079e-04 4.52833892e-04 2.97099958e-04 5.22922148e-04
 2.92262851e-04 2.20052284e-03 3.83532304e-04 6.32149739e-04
 2.65631266e-04 4.72368031e-04 7.89694271e-04 8.37232630e-04
 4.74232894e-04 5.64497691e-04 2.94122337e-03 2.63917959e-04
 5.04321725e-04 2.89069692e-04 6.90659396e-04 2.63585485e-04
 2.91107222e-04 2.83070273e-04 9.32592868e-04 5.66913333e-04
 4.08362849e-04 3.36558018e-03 2.76066020e-03 2.87497857e-04
 1.19049137e-03 4.52875210e-04 2.67604795e-04 4.02825302e-04
 7.05156479e-04 3.06836938e-04 2.78851619e-04 1.74883211e-03
 3.64067879e-04 6.07388648e-04 7.58277615e-04 7.57273054e-04
 5.29332460e-04 2.81986374e-04 2.82743794e-04 4.05569007e-04
 5.83592900e-04 7.04332174e-04 2.75563066e-04 4.31355869e-04
 4.42656246e-04 1.82385222e-03 5.22656398e-04 3.51338107e-04
 6.79176372e-04 1.03947054e-03 1.07116863e-03 1.48117322e-03
 2.65401165e-04 4.37486711e-04 1.09231708e-03 2.45871300e-04
 4.31703787e-04 3.31251709e-04 4.01178501e-04 1.48896621e-03
 3.56162782e+00 5.61537488e-04 4.83215948e-04 1.21070439e-03
 9.75200724e-04 4.01127609e-04 6.27011828e-04 2.66717361e-04
 2.95563009e-04 3.02506345e-04 1.44876067e-03 3.00789351e-04
 7.91476001e-04 4.39010114e-04 3.17796695e-04 2.98886372e-04
 3.91550120e-04 4.73051682e-04 6.31527738e-04 9.31808401e-04
 3.03824501e-04 3.79034197e-04 2.59611879e-04 6.49811471e-04
 2.87937183e-04 3.80029796e-03 5.59943324e-04 3.17157107e-04
 8.92633169e-04 4.14709560e-04 6.86773617e-01 2.48380968e-03
 2.90093477e-04 4.36361431e-04 2.59591922e-04 2.47050144e-04
 2.59685093e-04 3.01684047e-04 3.02507235e-04 5.18916384e-04
 2.84954521e-04 4.74586441e-04 3.04089487e-04 4.26933571e-04
 2.91491532e-03 3.17801973e-04 2.47929448e-04 9.06269529e-04
 3.94183033e-03 3.48312422e-04 5.50505354e-04 4.58721438e-04
 3.01822960e-03 2.52072586e-04 4.91210068e-04 5.05974483e-04
 4.22862009e-03 4.23452403e-04 1.21258864e-03 4.38637903e-04
 3.89214946e-04 3.75037377e-04 2.58273666e-04 2.52480169e-04
 2.55347542e-04 5.32865237e-04 1.62222644e-03 2.96881492e-04
 2.81489423e-04 5.10193227e-04 2.53749020e-04 2.53913870e-04
 7.48709098e-04 2.63204276e-04 3.58335313e-04 5.10050913e-04
 3.06070501e-04 2.69330637e-04 2.98364767e-04 3.28267579e-04
 3.54259807e-04 2.46176219e-04 3.72063464e-04 3.41298513e-04
 8.07934920e-04 3.66198609e-04 3.66939372e-03 3.67570572e-03
 3.90089084e-04 3.55853552e-04 2.60724438e-04 4.37018363e-04
 2.57131144e-04 5.06057744e-04 2.78888173e-04 4.97652487e-04
 3.39670036e-04 8.05463091e-04 3.51184750e-04 2.62145544e-04
 3.74375146e-04 3.82508222e-04 3.68841183e-03 4.19717652e-04
 4.55690629e-04 2.99726571e-04 2.89178308e-04 5.08179827e-04
 8.03202235e-04 2.94196485e-04 2.81584302e-04 9.80940709e-03
 6.16255126e-04 2.47946510e-04 2.63466143e-03 2.88074200e-04
 3.67563020e-04 3.03789671e-04 5.36379052e-04 4.62299803e-04
 2.80226710e-04 3.77007575e-04 7.47209568e-04 2.66495405e-04
 4.13772364e-04 7.21450302e-04 5.47056613e-04 7.37967987e-04
 4.27841466e-04 3.76690865e-04 4.37975867e-04 3.09631795e-03
 2.84687078e-04 2.55618180e-04 3.45388186e-03 2.93628277e-04
 2.61960005e-04 3.94595188e-04 5.89990377e-03 9.02572938e-04
 2.82159330e-04 4.91128582e-04 7.25198188e-04 3.00210510e-04
 3.85064819e-04 6.80836379e-04 6.44973092e-04 3.17077777e-04
 4.22945658e-04 4.01621666e-04 1.20344371e-03 1.07300818e-02
 2.66737814e-04 3.12277612e-04 3.20768810e-04 2.68265320e-04
 2.91625279e-03 2.79387537e-04 4.52872521e-04 3.01871767e-04
 6.26111818e-04 4.64855935e-04 3.62488941e-04 3.18872152e-04
 2.99944069e-04 4.51085526e-04 4.98276715e-04 3.30513367e-04
 3.31357117e-04 2.45649148e-04 3.97569437e-04 2.53923270e-04
 9.79669582e-04 2.97081472e-04 2.65974062e-04 2.83734175e-04
 2.45992907e-04 6.82390231e-04 3.25452929e-04 4.83621995e-04
 2.95153421e-04 2.82123615e-03 2.37722866e-03 2.60311484e-04
 2.78938163e-04 2.00983379e-03 9.78866811e-04 2.89492097e-04
 2.91580342e-04 3.87343121e-04 6.97254658e-04 4.36301857e-04
 3.17114784e-03 3.22186946e-04 3.85647265e-04 2.79761337e-03
 4.39186665e-04 4.94425448e-04 2.62787157e-04 5.30634894e-04
 2.94923616e-04 2.53532330e-04 7.41679657e-04 5.64260677e-04
 2.77355811e-04 7.28817918e-04 3.09293347e-04 2.98152988e-04
 6.19426368e-04 2.69393681e-04 2.93032510e-04 6.93512946e-04
 3.26959929e-04 4.69588837e-04 1.08769402e-03 5.98038443e-04
 2.76007963e-04 2.87637693e-04 4.54607469e-04 3.53033837e-04
 5.58269976e-03 3.65366914e-04 3.06447575e-03 2.58277589e-04
 8.63480494e-04 2.81210750e-04 6.35883420e-04 3.11540223e-04
 3.52162584e-04 3.06498378e-04 3.74048116e-04 2.54321920e-04
 3.80908175e-04 4.24951356e-03 2.84954354e-04 4.18691693e-01
 4.29198280e-04 2.34973565e-03 4.18652522e-04 4.86206706e-04
 1.26277968e-03 2.60487497e-04 1.23711073e-03 5.57395720e-04
 4.58085248e-04 3.92749045e-04 3.82798861e-04 2.85589944e-04
 2.96116397e-03 7.62436805e-04 3.62626773e-04 3.44899705e-04
 3.37202800e-04 2.92916378e-04 2.49218239e-04 7.95948605e-04
 3.12065059e-04 4.58312754e-04 4.16739905e-04 5.51838715e-04
 4.87625848e-04 4.50237635e-04 2.06161193e-03 3.32015810e-04
 4.83970721e-04 2.89200280e-04 3.67273098e-04 2.84789133e-04
 8.30385790e-04 2.72907722e-04 3.10285352e-04 9.78476629e-04
 5.03511525e-03]

In diesem Modell werden 963 Datensätze mit einer Cooks Distance > 4/n als mögliche Ausreiser identifiziert. Beispielsweise könnte versucht werden, den Einfluss dieser Datensätze auf das Modell durch Datentransformation (bspw. durch Verwendung des Logarithmus) zu verringern. Allerdings hilft diese Methode der Modellperformance genau so oft, wie sie ihr schadet. Daher wird hier auf eine Datentransformation verzichtet und die Datenpunkte stattdessen aus dem Trainings Set entfernt.

#drop outliers identified by Cook's distance
train_dataset2 = train_dataset.drop(train_dataset.index[out_d])
train_dataset2.info()
<class 'pandas.core.frame.DataFrame'>
Int64Index: 15419 entries, 14185 to 17879
Data columns (total 14 columns):
 #   Column                Non-Null Count  Dtype   
---  ------                --------------  -----   
 0   longitude             15419 non-null  float64 
 1   latitude              15419 non-null  float64 
 2   housing_median_age    15419 non-null  int64   
 3   total_rooms           15419 non-null  int64   
 4   total_bedrooms        15419 non-null  int64   
 5   population            15419 non-null  int64   
 6   households            15419 non-null  int64   
 7   median_income         15419 non-null  float64 
 8   median_house_value    15419 non-null  float64 
 9   ocean_proximity       15419 non-null  category
 10  people_per_household  15419 non-null  float64 
 11  rooms_per_household   15419 non-null  float64 
 12  bedrooms_per_people   15419 non-null  float64 
 13  bedrooms_per_rooms    15419 non-null  float64 
dtypes: category(1), float64(8), int64(5)
memory usage: 1.7 MB
fig = sns.scatterplot(data=train_dataset2, x='median_income', y='median_house_value')
_images/regression_99_0.png

Im Vergleich zum ursprünglichen Trainings Set, welches im Kapitel EDA ausführlich analysiert wurde, wurden einige X- und Y-Ausreiser entfernt. Dadurch hat sich die waagrechte Punktelinie bei 500000$ medin_house_value verkleinert und Ausreiser mit hohen X- und niedrigen Y-Werten wurden entfernt.

fig = sns.displot(data=train_dataset2, x='median_house_value')
_images/regression_101_0.png

Durch die Entfernung der Ausreiser hat sich das Histogramm von median_house_value etwas weiter normalisiert.

# Fit Model without outliers
lm = smf.ols(formula='median_house_value ~ median_income + ocean_proximity + bedrooms_per_rooms + people_per_household', data=train_dataset2).fit()
lm.summary()
OLS Regression Results
Dep. Variable: median_house_value R-squared: 0.739
Model: OLS Adj. R-squared: 0.739
Method: Least Squares F-statistic: 8725.
Date: Tue, 18 Jan 2022 Prob (F-statistic): 0.00
Time: 01:04:55 Log-Likelihood: -1.8952e+05
No. Observations: 15419 AIC: 3.790e+05
Df Residuals: 15413 BIC: 3.791e+05
Df Model: 5
Covariance Type: nonrobust
coef std err t P>|t| [0.025 0.975]
Intercept 3.341e+04 4404.814 7.585 0.000 2.48e+04 4.2e+04
ocean_proximity[T.INLAND] -6.1e+04 1098.946 -55.509 0.000 -6.32e+04 -5.88e+04
ocean_proximity[T.NEAR COAST] 6838.7281 1120.708 6.102 0.000 4642.008 9035.448
median_income 4.689e+04 365.489 128.285 0.000 4.62e+04 4.76e+04
bedrooms_per_rooms 3.58e+05 1.21e+04 29.558 0.000 3.34e+05 3.82e+05
people_per_household -2.492e+04 583.778 -42.686 0.000 -2.61e+04 -2.38e+04
Omnibus: 958.072 Durbin-Watson: 1.986
Prob(Omnibus): 0.000 Jarque-Bera (JB): 1206.262
Skew: 0.599 Prob(JB): 1.16e-262
Kurtosis: 3.664 Cond. No. 154.


Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

Durch die Entfernung der Ausreiser haben sich R-squared, adjusted R-squared und F-statistics deutlich verbessert.
AIC und BIC sind auch etwas kleiner geworden und haben sich somit verbessert.
Die p-values der Features bleiben unverändert.
Auch an den p-values für Jarque-Bera und Omnibus hat sich kaum etwas verändert.
Durbin-Watson liegt nun bei einem Wert <> 2. Generell ist bei einem Wert von > 2 eine negative Korrelation der Modellfehler anzunehmen.

fig = sm.graphics.plot_partregress_grid(lm)
fig.tight_layout(pad=0.2)
_images/regression_105_0.png

Die Partial Regression Plots visualisieren die Relation einer unabhängigen Variablen in Verbindung mit den anderen abhängigen Variablen im Modell in Bezug auf die abhängige Variable.
Die Plots zeigen eine negative Abhängigkeit von der Ausprägung I.INLAND von ocean_proximity und peope_per_household. Median_income und bedrooms_per_rooms scheinen eine positive Abhängigkeit auf median_house_value zu haben. Gemäß dem Diagramm hat die Variable ocean_proximity[T.NEAR COAST] keinen Einfluss auf median_income. Dies steht im Gegensatz zum p-value 0,00 dieser Variablen, nach dessen Interpretation die Variable einen Einfluss auf die abhängige Variable hat.

fig = sm.graphics.plot_regress_exog(lm, "ocean_proximity[T.NEAR COAST]")
fig.tight_layout(pad=1.0)
_images/regression_107_0.png

Auch diese Plots lassen darauf schließen, dass die Variable _ocean_proximity[T.NEAR COAST] das Modell nicht verbessert. Sowohl Partial Regression Plot als auch CCPR Plot visualisieren eine Nicht-Linearität zwischen der Variablen und median_house_value.
Die Residuals scheinen dagegen gleichmäßig verteilt zu sein.

fig = sm.graphics.plot_regress_exog(lm, "ocean_proximity[T.INLAND]")
fig.tight_layout(pad=1.0)
_images/regression_109_0.png

Die negative Korrelation ist sowohl im Partial Regression Plot als auch im CCPR Plot zu erkennen. Die Residuals sind gleichmäßig verteilt und lassen daher auf Homoskedastizität rückschließen. Das Diagramm zu Y and Fitted vs. X ist aufgrund der vielen Datenpunkte schwer zu analysieren.

fig = sm.graphics.plot_regress_exog(lm, "people_per_household")
fig.tight_layout(pad=1.0)
_images/regression_111_0.png

Auch hier ist eine negative Abhängigkeit der abhängigen zur unabhängigen Variablen zu erkennen. Linearität scheint also vorzuliegen. Allerdings ist im Residual Plot keine konstante Varianz der Fehlerstreuung zu erkennen. Mit steigendem X-Wert verringert sich die Fehlerstreuung und liegt nur noch im positiven Bereich. Mit zunehmender people_per_household wird tendenziell ein höherer median_house_value mit dem Modell hervorgesagt, als tatsächlich vorliegt. Dies spiegelt sich auch im Plot zu Y and Fitted vs. X wieder. Bei dieser Variablen liegt somit Heteroskedastizität vor. Die Ursache hierfür könnte darin liegen, dass im höheren Wertebereich von people_per_household deutlich weniger Datensätze vorliegen als im niedrigeren Bereich und das Modell somit nicht ausreichend für diese Wertebereiche trainiert wurde. Womöglich bildet auch eine polynomiale Modellierung für diese Variable die Relation der Variablen besser ab.

fig = sm.graphics.plot_regress_exog(lm, "median_income")
fig.tight_layout(pad=1.0)
_images/regression_113_0.png

Für median_income ist die Linearität im Partial Regression Plot und CCPR Plot deutlich zu erkennen. Jedoch nimmt auch hier die Streuung der Residuals im höheren X-Wertebereich ab.

# write out a dict with the mapping of old to new
remap_cat_dict = {
    'NEAR COAST': 'NEAR COAST',
    'INLAND': 'INLAND',
    '<1H OCEAN': 'NEAR COAST' }

train_dataset.ocean_proximity = train_dataset.ocean_proximity.map(remap_cat_dict).astype('category')

Model 2 - OLS Regression

Nachfolgend wird ein weiteres Modell mit zusätzlicher Verwendung der Variablen housing_median_age trainiert.

# Fit Model
lm2 = smf.ols(formula='median_house_value ~ median_income + ocean_proximity + bedrooms_per_rooms + people_per_household + housing_median_age', data=train_dataset).fit()
lm2.summary()
OLS Regression Results
Dep. Variable: median_house_value R-squared: 0.627
Model: OLS Adj. R-squared: 0.627
Method: Least Squares F-statistic: 5486.
Date: Tue, 18 Jan 2022 Prob (F-statistic): 0.00
Time: 01:08:50 Log-Likelihood: -2.0563e+05
No. Observations: 16340 AIC: 4.113e+05
Df Residuals: 16334 BIC: 4.113e+05
Df Model: 5
Covariance Type: nonrobust
coef std err t P>|t| [0.025 0.975]
Intercept -7.016e+04 4213.120 -16.652 0.000 -7.84e+04 -6.19e+04
ocean_proximity[T.NEAR COAST] 6.563e+04 1348.977 48.655 0.000 6.3e+04 6.83e+04
median_income 4.445e+04 406.063 109.460 0.000 4.37e+04 4.52e+04
bedrooms_per_rooms 2.998e+05 1.29e+04 23.243 0.000 2.74e+05 3.25e+05
people_per_household -1.096e+04 428.532 -25.570 0.000 -1.18e+04 -1.01e+04
housing_median_age 1005.0830 46.230 21.741 0.000 914.467 1095.699
Omnibus: 4984.528 Durbin-Watson: 1.968
Prob(Omnibus): 0.000 Jarque-Bera (JB): 27788.201
Skew: 1.356 Prob(JB): 0.00
Kurtosis: 8.785 Cond. No. 766.


Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

R-squared und adj. R-squared haben sich im Vergleich zum ersten Modell verbessert.
Auch F-statistics hat sich durch die Entfernung der Ausprägung von ocean_proximity und Aufnahme des Features median_housing_age verbessert.
Alle übrigen Kennzahlen haben sich kaum verändert.

# obtain Cook's distance 
lm_cooksd = lm2.get_influence().cooks_distance[0]

# get length of df to obtain n
n = len(train_dataset["median_income"])
print('n = ', n)

# calculate critical d
critical_d = 4/n
print('Critical Cooks distance:', critical_d)

# identification of potential outliers with leverage
out_d = lm_cooksd > critical_d

# output potential outliers with leverage
print(train_dataset.index[out_d], "\n", 
    lm_cooksd[out_d])
n =  16340
Critical Cooks distance: 0.0002447980416156671
Int64Index([15778,  6776, 10005, 11511,  7263, 15723, 14756, 15687,  4638,
            14806,
            ...
            15809,  8866,   199, 18498,    59, 10749,  8857, 13154,  5470,
            16642],
           dtype='int64', length=834) 
 [1.34620956e-03 3.28749861e-04 2.57215207e-04 3.35039575e-04
 6.51518260e-04 3.61919469e-04 7.91550903e-03 6.11015995e-04
 6.49112668e-04 4.17870132e-04 8.15554753e-04 1.30731202e-03
 1.20416904e-03 3.73036827e-04 5.62180943e-04 7.70163204e-04
 5.46624084e-04 2.53753610e-04 3.01304709e-04 5.38305099e-04
 3.35811574e-04 2.54781398e-04 4.29838436e-04 3.35159404e-04
 3.69843859e-04 3.46408581e-04 2.57119718e-04 5.47776435e-02
 4.46822263e-04 4.27225287e-04 3.34569001e-04 8.73800145e-04
 5.12264946e-04 1.65774076e-03 2.06412271e-03 1.04840879e-03
 2.72011719e-04 5.81620070e-04 5.20076218e-04 5.13152695e-03
 5.17172666e-04 9.28482089e-04 4.71110244e-03 3.95928870e-04
 2.95769203e-04 2.62218028e-04 9.37446313e-04 2.83036311e-04
 3.06475831e-04 4.26676622e-04 3.01809213e-04 5.76534941e-04
 3.38711782e-04 3.70526120e-04 4.21168257e-04 1.19851082e-03
 3.09081908e-04 6.12841623e-04 3.59409818e-04 3.20014681e-04
 3.94920532e-04 9.53529426e-04 2.57227742e-04 8.06320199e-04
 2.55183702e-04 4.39410335e-04 6.62378271e-04 5.62770048e-04
 3.31933928e-04 7.31914431e-04 1.00363583e-03 5.29889822e-04
 4.33576686e-04 4.95675084e-04 8.89908119e-04 2.83041383e-04
 8.55849228e-04 2.16001576e-03 3.09497542e-04 6.81234101e-04
 7.08301229e-04 6.74865681e-04 3.14174140e-04 9.77863893e-04
 2.98076323e-04 8.43691811e-04 1.18174422e-03 9.89622636e-04
 5.13546336e-04 9.31650342e-04 1.20073880e-03 3.56089239e-04
 6.11493487e-04 4.13533793e-04 1.64662507e-03 3.81249929e-04
 9.84393006e-04 6.50545586e-04 2.18054214e-03 2.78936483e-04
 2.52932373e-04 2.52777940e-04 2.66069166e-04 5.83492758e-04
 3.39562190e-04 6.77350435e-04 6.48707359e-04 4.21307444e-04
 2.91646476e-04 2.47710633e-04 2.53301281e-04 3.56535739e-04
 7.43918045e-04 3.77300968e-04 8.44956605e-04 1.30709420e-03
 2.85992118e-04 4.39369242e-04 3.80927914e-04 3.82280250e-04
 4.09153132e-04 3.44268183e-04 5.46423976e-04 4.20173053e-04
 4.28627490e-02 7.35095993e-04 4.24027988e-04 7.22312913e-04
 1.32786505e-03 6.99537001e-04 2.50626133e-04 5.46342106e-04
 4.94243518e-04 3.65875470e-04 3.13187264e-04 2.70271254e-04
 3.25703718e-04 6.93446296e-04 2.65491825e-04 3.49024635e-04
 4.08335882e-04 4.05092100e-04 2.45889472e-04 4.39110074e-04
 2.85502846e-04 9.38717658e-04 1.34362145e-03 4.80589164e-03
 3.98106969e-04 3.22240865e-04 3.15091984e-04 2.56911300e-04
 3.46225434e-04 4.22497484e-04 3.46518228e-04 6.70520360e-04
 1.22863087e-03 7.43148088e-04 2.54542565e-04 2.56665208e-04
 1.40384851e-03 3.28558370e-03 3.11047780e-04 3.43185408e-04
 1.43474503e-03 7.15000996e-04 4.19995923e-03 7.80899608e-04
 3.75125809e-04 2.79783608e-04 3.39099740e-04 1.19937066e-03
 4.88403117e-03 3.65830328e-04 5.86482204e-04 2.56515079e-04
 4.23646462e-03 2.54733684e-04 4.59159574e-04 8.87961889e-03
 3.74045810e-03 4.00718231e-03 5.46171186e-04 1.50078936e-03
 3.56173603e-04 4.34240443e-04 2.61806316e-04 3.31107538e-04
 1.18033323e-03 2.47950011e-04 3.32658507e-04 3.22346407e-04
 2.61654415e-04 6.16450930e-04 9.91135344e-04 1.06921629e+01
 2.82232779e-04 3.13847258e-04 3.09257186e-04 3.06408629e-04
 2.56403994e-04 1.05400214e-03 4.00527883e-02 2.84546161e-04
 2.53531661e-04 2.11139733e-03 3.83398773e-04 4.32087727e-04
 1.45600603e-03 6.41920570e-04 6.14620185e-04 3.14498312e-04
 2.97336263e-04 3.07533280e-02 8.85053202e-04 4.54537558e-04
 4.20079773e-04 3.73657618e-04 2.85897551e-04 1.45010004e-03
 4.47302087e-04 2.86417016e-04 3.86170479e-04 3.64811817e-04
 9.04250746e-04 8.96897085e-04 3.31137357e-03 4.58857362e-04
 6.13777854e-04 2.81677164e-04 3.09589735e-04 5.65060691e-04
 9.58478200e-04 3.69278949e-04 1.84718845e-03 2.58189555e-04
 3.25494577e-04 4.15475232e-04 4.83547910e-03 2.45981050e-04
 2.75498461e-04 7.50737541e-04 4.33679600e-03 3.02966513e-04
 1.16482267e-03 4.43903600e-04 4.39257693e-04 5.23686883e-04
 5.43692280e-04 6.87454121e-04 3.78535731e-04 2.49192273e-04
 4.34747272e-04 3.39985468e-04 3.30260839e-04 2.66248512e-03
 3.09325047e-04 3.81360279e-04 6.38501803e-04 7.93860341e-04
 3.17777251e-03 4.86447450e-03 4.16377804e-04 3.53579296e-04
 2.94406173e-04 2.61202267e-04 3.66836210e-04 2.62061270e-04
 3.38031251e-04 4.27777058e-04 9.13378989e-04 3.41456845e-04
 3.64618565e-03 4.72926716e-04 1.77649059e-03 5.02697205e-04
 3.67504825e-03 1.71681800e-03 5.55103694e-04 2.91051755e-04
 3.55591228e-04 4.59132619e-04 2.07003198e-03 3.97897537e-04
 3.05201888e-04 6.45379090e-04 3.27399025e-04 4.93574804e-04
 4.46207746e-03 1.06687385e-03 1.27664267e-03 5.98190815e-04
 6.69362545e-03 6.47563661e-04 1.41449488e-03 3.35907072e-04
 4.31922192e-04 9.62879839e-04 5.99119839e-04 5.94298294e-04
 2.72450549e-04 1.33284800e-03 2.22449093e-03 2.94392265e-04
 3.60548234e-04 2.46831910e-04 5.21180548e-03 3.74112820e-04
 3.00529729e-04 1.83317278e-03 2.10050482e-03 2.46498466e-04
 5.52625587e-04 5.69750888e-04 2.70235120e-04 4.73789257e-04
 6.37396918e-04 1.33110610e-03 3.97863971e-04 3.91811176e-04
 5.56690305e-04 3.49347749e-04 3.60971784e-04 2.55520221e-04
 3.55568409e-04 3.07803326e-04 1.03831598e-03 2.82099059e-04
 3.35719194e-04 2.63385336e-04 7.30688902e-04 2.46434487e-04
 3.32206636e-04 5.71495033e-04 4.55220555e-04 4.58525682e-04
 2.51076913e-04 7.81748699e-03 6.88520513e-04 4.47783205e-04
 1.38917952e-03 4.58028391e-04 8.53315502e-04 7.98043674e-04
 3.06666704e-04 3.57944205e-04 2.97907283e-04 6.60775311e-04
 3.15356513e-04 2.64416584e-04 4.01810556e-04 4.11109209e-04
 1.17853473e-03 3.03560371e-04 7.12310609e-04 7.81473143e-04
 2.45443831e-04 4.31752702e-04 3.55261927e-04 4.28943588e-04
 1.18408870e-03 1.26262767e-03 3.30624525e-04 3.11377904e-04
 7.37247025e-04 4.34417931e-04 3.13091517e-04 3.86335627e-04
 3.27376743e-04 1.69310467e-03 3.75865197e-04 2.94169921e-04
 4.09415944e-03 5.42417148e-04 5.02146968e-03 2.72579606e-04
 3.40573948e-04 1.36489576e-03 3.93579498e-03 5.59166602e-04
 3.11084962e-04 3.49829003e-04 3.02065944e-03 4.74636086e-04
 8.17112634e-04 5.40381997e-04 2.21043242e-03 2.99376796e-04
 4.78605704e-03 2.65136147e-04 3.62158887e-04 5.37258826e-04
 7.16663410e-04 3.34015704e-04 2.37806674e-03 3.87874124e-04
 4.60018714e-03 4.77441784e-04 3.16011510e-04 8.63052283e-04
 6.69780131e-04 2.87389645e-04 4.19416247e-04 5.04048059e-04
 2.01334240e-03 3.63296782e-04 3.86119046e-03 6.21399855e-04
 6.91604511e-04 3.00672616e-04 3.07340296e-04 8.51472390e-04
 1.39838656e-03 6.99686628e-04 3.88329899e-04 1.44844319e-01
 4.06391695e-03 3.81700667e-04 4.36672272e-04 6.56546831e-04
 5.04080380e-04 1.40540348e-03 4.24763158e-04 5.29649070e-03
 2.49851273e-04 1.37361535e-03 5.45285648e-04 7.92267019e-04
 8.89794507e-04 2.05090829e-03 3.00435326e-04 3.57210150e-04
 4.07742500e-04 1.18352972e-03 2.43728565e-03 2.82272469e-04
 3.89905770e-04 2.25566831e-03 3.03114176e-04 1.30068085e-02
 2.58994479e-04 5.53658856e-04 6.12830952e-04 5.01723767e-04
 8.12629182e-04 3.77644531e-04 1.55498141e-03 3.05627947e-04
 5.94422552e-03 3.57403852e-04 8.78405940e-03 3.13205038e-04
 4.27447728e-04 7.44573473e-04 5.15549175e-03 7.62643519e-03
 3.32130683e-04 2.61203844e-04 1.29768456e-03 2.48558136e-04
 7.25434452e-04 3.67663029e-03 4.20029578e-04 4.90839693e-04
 2.73133315e-04 3.38691084e-04 1.89005180e-03 4.28670600e-04
 2.55054527e-04 3.18588498e-04 7.39903481e-04 2.31834749e-03
 3.67185876e-04 3.65229676e-04 5.75161805e-04 3.15816880e-04
 4.74515715e-04 5.13067456e-04 6.89572009e-04 2.77145250e-04
 4.70737869e-04 2.64762172e-04 5.27422472e-04 9.84707517e-04
 9.43212824e-03 1.70314323e-03 4.62703809e-04 3.92885812e-04
 4.20675602e-04 1.95409614e-03 2.88329437e-04 1.32460983e-03
 6.43856512e-04 4.88842954e-04 2.69788415e-04 3.59022621e-03
 3.48583713e-04 1.65698848e-03 2.23761619e-03 3.61044126e-04
 7.20461757e-04 3.31346884e-04 2.44398609e-03 7.30058445e-04
 4.79412548e-04 6.82860694e-04 8.29685275e-04 1.15423848e-03
 4.92749800e-03 8.76630693e-04 2.63359256e-04 5.61268507e-04
 7.13278888e-04 1.70939726e-03 2.71500955e-04 8.65533017e-04
 8.66559196e-04 6.66489288e-04 3.82956960e-04 3.70377201e-04
 3.62988149e-04 2.48974052e-04 2.90398752e-04 3.77658563e-04
 2.49335881e-03 2.53154039e-04 2.66962782e-04 9.60055604e-04
 5.36687000e-04 6.75338838e-04 2.07363805e-03 7.50365981e-04
 3.87624234e-04 4.74190135e-03 2.87536148e-04 3.54576216e-04
 8.61051300e-04 3.07458358e-04 3.37358573e-04 3.19957211e-04
 2.92073680e-04 6.06543953e-04 5.92986873e-04 3.08890432e-04
 3.43232872e-03 3.15607985e-03 3.87409151e-04 1.20192208e-03
 4.33069594e-04 7.01263187e-04 3.66029598e-04 1.51065092e-03
 5.10965707e-04 9.46617435e-04 1.33914359e-03 5.42132639e-04
 2.55992635e-04 2.79273013e-04 4.33360234e-04 9.02970286e-04
 6.22861701e-04 5.34084933e-04 4.66705558e-04 3.30393525e-03
 2.95158501e-04 5.58204627e-04 2.76748391e-04 9.28388580e-04
 9.60918152e-04 8.80922972e-04 1.07105882e-03 3.69483391e-04
 4.64798647e-04 2.73824841e-04 9.66985466e-04 2.83190149e-04
 5.82375269e-04 2.60887897e-04 3.25360840e-04 1.15909394e-03
 4.02187581e+00 5.43475220e-04 4.02644115e-04 8.64339222e-04
 8.18907455e-04 3.24992180e-04 7.38184531e-04 2.51052121e-04
 6.42272514e-04 3.04544105e-04 3.63665615e-04 1.02912398e-03
 2.87378521e-04 7.45509328e-04 3.26415591e-04 2.52550349e-04
 3.40003285e-04 4.24117891e-04 4.90316626e-04 4.86412402e-04
 6.40754977e-04 3.39882911e-04 6.81558602e-04 3.32647839e-04
 3.38165023e-03 4.58021592e-04 6.59793292e-04 5.37127516e-04
 7.46891438e-01 2.20109675e-03 3.91297892e-04 4.89698239e-04
 2.60150085e-04 3.59658963e-04 3.11584443e-04 3.02795328e-04
 4.18784806e-04 2.78422616e-04 4.22763551e-04 2.65366015e-04
 2.65452108e-04 2.68691754e-04 2.98144020e-04 3.86472350e-03
 3.38594356e-04 8.99398215e-04 5.10549736e-03 4.16733375e-04
 4.03199774e-04 7.18913911e-04 4.79645149e-04 5.40162714e-03
 3.87982015e-04 6.20419951e-04 4.25703669e-03 4.71147289e-04
 1.04261716e-03 3.70599935e-04 4.27817891e-04 6.19133596e-04
 2.53943173e-04 3.41160201e-04 1.66679303e-03 6.29287626e-04
 5.17842252e-04 3.12862593e-04 7.27747807e-04 2.61667141e-04
 4.64546973e-04 4.43050694e-04 2.89239694e-04 3.40418715e-04
 6.36218958e-04 2.98524754e-04 8.09693242e-04 1.15337376e-03
 2.63411919e-04 3.58669103e-03 2.98727727e-03 2.53558984e-04
 3.42658703e-04 5.46202372e-04 3.09921088e-04 4.24199927e-04
 2.66680028e-04 2.52166862e-04 4.51032150e-04 3.98717184e-04
 7.78799300e-04 4.43084627e-04 2.66797596e-04 2.70016395e-04
 3.52045620e-04 6.71406313e-04 3.51355939e-03 4.71552792e-04
 2.67398242e-04 3.26641206e-04 4.89193724e-04 5.54156290e-04
 3.55441606e-04 8.70449188e-03 1.01979909e-03 4.47454853e-03
 2.67577330e-04 6.03933837e-04 4.19977124e-04 3.30883814e-04
 2.59675243e-04 7.82055746e-04 5.70203443e-04 2.82550451e-04
 7.75940421e-04 5.66976263e-04 7.23078868e-04 5.58070994e-04
 3.80964154e-04 4.70836891e-04 4.63370304e-03 2.93061722e-04
 2.82097527e-04 5.25394787e-04 4.39058608e-03 3.02814301e-04
 3.93885720e-04 2.72975994e-04 6.30941339e-04 6.34256527e-03
 1.93264396e-03 2.63895051e-04 3.90119304e-04 6.49958891e-04
 3.39334191e-04 2.80016960e-04 6.70981971e-04 7.96887450e-04
 2.61031562e-04 2.99531345e-04 3.92412896e-04 9.15393354e-04
 9.65937206e-03 2.59496858e-04 3.04018109e-04 3.35152782e-04
 4.47545315e-03 3.28849297e-04 3.09466337e-04 2.57128750e-04
 4.79134179e-04 2.84444223e-04 7.04680802e-04 3.08216566e-04
 3.19066173e-04 3.31656683e-04 4.99721752e-04 3.88382147e-04
 4.62130992e-04 3.38338725e-04 3.33453243e-04 3.36783921e-04
 3.53353487e-04 4.69560563e-04 1.01286041e-03 2.65157668e-04
 4.44166808e-04 5.14610488e-04 3.25086702e-04 4.05530866e-03
 2.61966032e-03 3.81661774e-04 1.26797417e-03 1.00604063e-03
 2.47097398e-04 2.54264288e-04 3.76542325e-04 5.62680829e-04
 4.39275303e-04 3.78306069e-03 4.09577971e-04 4.23492354e-03
 4.67446416e-04 4.53267869e-04 2.61729406e-04 4.40048926e-04
 7.43015952e-04 5.68428522e-04 2.68099048e-04 3.75476250e-04
 6.99763254e-04 3.25127465e-04 2.44953335e-04 5.51778636e-04
 3.10795171e-04 3.74624073e-04 6.56877859e-04 4.18323970e-04
 4.79464330e-04 7.98927851e-04 3.65242322e-04 2.69035991e-04
 3.39453306e-04 4.82734528e-04 3.59781628e-04 6.59773820e-03
 3.12038525e-04 2.57552160e-04 3.11064867e-03 5.09224651e-04
 7.60255237e-04 2.55022118e-04 5.33067322e-04 3.11642549e-04
 2.90908224e-04 6.34810463e-04 4.38153846e-04 5.71802541e-03
 6.88659684e-04 3.99323486e-01 5.10145404e-04 2.98228342e-03
 3.12344844e-04 5.30124261e-04 9.11032421e-04 1.47777500e-03
 1.02556616e-03 4.81595818e-04 4.69679318e-04 2.80467578e-04
 2.75736314e-03 1.01866312e-03 4.80632918e-04 3.61741717e-04
 2.65964061e-04 2.66946243e-04 1.07338565e-03 2.64243381e-04
 3.62524321e-04 2.64787790e-04 3.70798360e-04 2.91009002e-04
 6.70681350e-04 3.49234458e-04 3.54200151e-04 1.38738166e-03
 3.21375424e-04 3.65064441e-04 3.17517938e-04 4.03136086e-04
 2.73400173e-04 8.99485876e-04 2.73811960e-04 3.31912000e-04
 9.43051451e-04 5.70143496e-03]

Auch in Modell lm2 werden mögliche outliers und Punkte mit high leverage identifiziert und aus dem Datenset entfernt.

#drop outliers identified by Cook's distance
train_dataset2 = train_dataset.drop(train_dataset.index[out_d])
train_dataset2.info()
<class 'pandas.core.frame.DataFrame'>
Int64Index: 15506 entries, 14185 to 17879
Data columns (total 14 columns):
 #   Column                Non-Null Count  Dtype   
---  ------                --------------  -----   
 0   longitude             15506 non-null  float64 
 1   latitude              15506 non-null  float64 
 2   housing_median_age    15506 non-null  int64   
 3   total_rooms           15506 non-null  int64   
 4   total_bedrooms        15506 non-null  int64   
 5   population            15506 non-null  int64   
 6   households            15506 non-null  int64   
 7   median_income         15506 non-null  float64 
 8   median_house_value    15506 non-null  float64 
 9   ocean_proximity       15506 non-null  category
 10  people_per_household  15506 non-null  float64 
 11  rooms_per_household   15506 non-null  float64 
 12  bedrooms_per_people   15506 non-null  float64 
 13  bedrooms_per_rooms    15506 non-null  float64 
dtypes: category(1), float64(8), int64(5)
memory usage: 1.7 MB
# Fit Model
lm2 = smf.ols(formula='median_house_value ~ median_income + ocean_proximity + bedrooms_per_rooms + people_per_household + housing_median_age', data=train_dataset2).fit()
lm2.summary()
OLS Regression Results
Dep. Variable: median_house_value R-squared: 0.743
Model: OLS Adj. R-squared: 0.743
Method: Least Squares F-statistic: 8956.
Date: Tue, 18 Jan 2022 Prob (F-statistic): 0.00
Time: 01:09:17 Log-Likelihood: -1.9073e+05
No. Observations: 15506 AIC: 3.815e+05
Df Residuals: 15500 BIC: 3.815e+05
Df Model: 5
Covariance Type: nonrobust
coef std err t P>|t| [0.025 0.975]
Intercept -5.584e+04 4076.067 -13.700 0.000 -6.38e+04 -4.79e+04
ocean_proximity[T.NEAR COAST] 5.731e+04 1067.209 53.703 0.000 5.52e+04 5.94e+04
median_income 4.88e+04 370.848 131.580 0.000 4.81e+04 4.95e+04
bedrooms_per_rooms 3.69e+05 1.19e+04 30.882 0.000 3.46e+05 3.92e+05
people_per_household -2.592e+04 562.240 -46.097 0.000 -2.7e+04 -2.48e+04
housing_median_age 924.7246 36.788 25.137 0.000 852.616 996.833
Omnibus: 1156.804 Durbin-Watson: 1.988
Prob(Omnibus): 0.000 Jarque-Bera (JB): 1572.169
Skew: 0.645 Prob(JB): 0.00
Kurtosis: 3.876 Cond. No. 911.


Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

Die Modellperformance liegt nun bei 74,3% und ist somit bisher der beste erzielte Wert der Kennzahl R-squared.
F-statistics hat sich ebenfalls verbessert, auch im Vergleich zu Model 1.
Alle anderen Werte liegen in einem ähnlichen Wertebereich.

fig = sm.graphics.plot_partregress_grid(lm2)
fig.tight_layout(pad=0.2)
_images/regression_125_0.png
fig = sm.graphics.plot_regress_exog(lm2, "housing_median_age")
fig.tight_layout(pad=1.0)
_images/regression_126_0.png

Das Partial Regression Plot und CCPR Plot zeigen eine geringe lineare Korrelation von housing_median_age in Verbindung mit den anderen unabhängigen Variablen auf median_house_value.
Das Residual Plot zeigt eine deutliche Homoskedastizität.

fig = sm.graphics.plot_regress_exog(lm2, "ocean_proximity[T.NEAR COAST]")
fig.tight_layout(pad=1.0)
_images/regression_128_0.png

Die Plots von ocean_proximity[T.NEAR COAST] zeigen keine große Veränderung zu Modell lm.

In der linearen Regression muss eine Colinearität der unabhängigen Variablen unbedingt vermieden werden. Mithilfe des Variance Inflation Factors (VIF) kann der Grad der Colinearität der im Modell verwendeten Variablen berechnet werden. Der kleinstmögliche Wert ist 1. Problematisch wird es ab Werten > 5.

from statsmodels.stats.outliers_influence import variance_inflation_factor
from statsmodels.tools.tools import add_constant

# choose features and add constant
features = add_constant(train_dataset2[['median_income', 'bedrooms_per_rooms', 'people_per_household', 'housing_median_age']])
# create empty DataFrame
vif = pd.DataFrame()
# calculate vif
vif["VIF Factor"] = [variance_inflation_factor(features.values, i) for i in range(features.shape[1])]
# add feature names
vif["Feature"] = features.columns

vif.round(2)
C:\ProgramData\Anaconda3\lib\site-packages\statsmodels\tsa\tsatools.py:142: FutureWarning:

In a future version of pandas all arguments of concat except for the argument 'objs' will be keyword-only
VIF Factor Feature
0 82.06 const
1 1.80 median_income
2 1.79 bedrooms_per_rooms
3 1.01 people_per_household
4 1.03 housing_median_age

Der Wert der Konstanten kann ignoriert werden. Ansonsten ist die Colinearität der Variablen eher gering.

Die Modelle lm und lm2 sind sehr ähnlich. R-squared ist bei lm2 ca. 1% besser, dafür hat sich der Wert von F-statistics von lm zu lm2 um ca. 1000 verschlechtert. AIC und BIC haben sich nicht wesentlich verändert. Die Werte sind von lm zu lm2 sogar geringfügig kleiner (also besser) geworden, obwohl für zusätzliche Features ein Penalty Term einberechnet wird, um overfitting zu vermeiden. Das Partial Regression Plot der Variablen median_housing_age zeigt zwar nur eine geringe Linearität, dennoch scheint die Aufnahme der Variablen das Modell etwas zu verbessern.

Model 2 - Lasso Regression

Nachfolgend wird eine Lasso Regression mit statsmodels fit_regularized() Methode trainiert. Der Wert von alpha wird auf 5 gesetzt.

lasso1 = smf.ols(formula='median_house_value ~ median_income + ocean_proximity + bedrooms_per_rooms + people_per_household + housing_median_age', data=train_dataset2).fit_regularized(method='sqrt_lasso', alpha=5.0, L1_wt=1.0, max_iter=1000)
print("Parameters: ", lasso1.params)
Parameters:  Intercept                        -50616.495719
ocean_proximity[T.NEAR COAST]     57809.976303
median_income                     48394.058880
bedrooms_per_rooms               351976.704469
people_per_household             -26004.829908
housing_median_age                  918.971072
dtype: float64

Da die Summary() Methode für statsmodels fit_regularized() Medelle nicht verfügbar ist, müssen wichtige Kennzahlen selbst berechnet werden.

# calculate the mean 
train_dataset2['average'] = train_dataset2["median_house_value"].mean()

# calculate error (observation - average) and assign it to dataframe
train_dataset2 = train_dataset2.assign(error = (train_dataset2['median_house_value'] - train_dataset2['average']))

# calculate squared error and assign it to dataframe
train_dataset2 = train_dataset2.assign(error_sq = (train_dataset2['median_house_value'] - train_dataset2['average'])**2)


# obtain predicted value from model lasso
train_dataset2['pred'] = lasso1.predict()

# obtain the residuals from statsmodel (resid)
train_dataset2['error_2'] = train_dataset2['median_house_value'] - train_dataset2['pred']

# square the residuals 
train_dataset2['error_sq_2'] = train_dataset2['error_2']**2

# show df
train_dataset2.head(5)
longitude latitude housing_median_age total_rooms total_bedrooms population households median_income median_house_value ocean_proximity people_per_household rooms_per_household bedrooms_per_people bedrooms_per_rooms average error error_sq pred error_2 error_sq_2
14185 -117.08 32.70 37 2176 418 1301 375 2.8750 98900.0 NEAR COAST 3.469333 5.802667 0.321291 0.192096 196674.726944 -97774.726944 9.559897e+09 157722.078408 -58822.078408 3.460037e+09
6125 -117.91 34.11 20 3158 684 2396 713 3.5250 153000.0 NEAR COAST 3.360449 4.429173 0.285476 0.216593 196674.726944 -43674.726944 1.907482e+09 185009.672924 -32009.672924 1.024619e+09
14095 -117.10 32.75 11 2393 726 1905 711 1.3448 91300.0 NEAR COAST 2.679325 3.365682 0.381102 0.303385 196674.726944 -105374.726944 1.110383e+10 119491.512237 -28191.512237 7.947614e+08
14359 -117.22 32.74 52 1260 202 555 209 7.2758 345200.0 NEAR COAST 2.655502 6.028708 0.363964 0.160317 196674.726944 148525.273056 2.205976e+10 394457.593271 -49257.593271 2.426310e+09
18004 -121.99 37.29 32 2930 481 1336 481 6.4631 344100.0 NEAR COAST 2.777547 6.091476 0.360030 0.164164 196674.726944 147425.273056 2.173421e+10 334928.406586 9171.593414 8.411813e+07
# calculate sum of squared error (which is in case of the mean the total error)
TSS = train_dataset2.error_sq.sum()
# print output
print('Sum of squared error (TSS) of model 1:', TSS)

# Sum of squared residuals (SS_R)
SSR = train_dataset2['error_sq_2'].sum()
print('Sum of squared residuals (SSR) of model 1:',SSR)
Sum of squared error (TSS) of model 1: 170542426261911.88
Sum of squared residuals (SSR) of model 1: 43859574696952.75
# Explained sum of squares  (SS_M = TSS - SS_R)
SSM = TSS - SSR
print('SSM of Lasso Model: ', SSM)

# R_Squared: explained sum of squared residuals
R_squared = SSM / TSS
print('R squared Lasso Model:', R_squared)
SSM of Lasso Model:  126682851564959.12
R squared Lasso Model: 0.7428230871443388

R-squared ist durch die Regularisierung nochmal etwas besser geworden.
Für einen optimalen Wert von alpha müsste eigentlich eine Cross-Validation durchgeführt werden. Darauf wird an dieser Stelle jedoch verzichtet.

Model 3 - Cubic Spline

Nun wird eine Cubic Spline unter Verwendung von Patsy und Statsmodels trainiert. Eine Verwendung von mehreren unabhängigen Variablen ist dabei nicht ohne Weiteres möglich. Daher wird nachfolgend lediglich die Variable mit der höchsten Korrelation zu median_house_value verwendet. Dabei handelt es sich um median_income.

X_train = train_dataset2['median_income']
y_train = train_dataset2['median_house_value']
X_test = test_dataset['median_income']
y_test = test_dataset['median_house_value']
from patsy import dmatrix

# Generating cubic spline with 3 knots at 25, 40 and 60
transformed_x = dmatrix(
            "bs(train, knots=(3,6,9), degree=3, include_intercept=False)", 
                {"train": X_train},return_type='dataframe')
# Fitting generalised linear model on transformed dataset
spline1 = sm.GLM(y_train, transformed_x).fit()
spline1.predict()
spline1.summary()
Generalized Linear Model Regression Results
Dep. Variable: median_house_value No. Observations: 15506
Model: GLM Df Residuals: 15499
Model Family: Gaussian Df Model: 6
Link Function: identity Scale: 4.6148e+09
Method: IRLS Log-Likelihood: -1.9452e+05
Date: Tue, 18 Jan 2022 Deviance: 7.1525e+13
Time: 00:43:42 Pearson chi2: 7.15e+13
No. Iterations: 3
Covariance Type: nonrobust
coef std err z P>|z| [0.025 0.975]
Intercept 1.081e+05 8237.478 13.122 0.000 9.19e+04 1.24e+05
bs(train, knots=(3, 6, 9), degree=3, include_intercept=False)[0] -3.748e+04 1.09e+04 -3.448 0.001 -5.88e+04 -1.62e+04
bs(train, knots=(3, 6, 9), degree=3, include_intercept=False)[1] 7.416e+04 7398.863 10.023 0.000 5.97e+04 8.87e+04
bs(train, knots=(3, 6, 9), degree=3, include_intercept=False)[2] 1.575e+05 9980.819 15.782 0.000 1.38e+05 1.77e+05
bs(train, knots=(3, 6, 9), degree=3, include_intercept=False)[3] 4.47e+05 1.06e+04 41.988 0.000 4.26e+05 4.68e+05
bs(train, knots=(3, 6, 9), degree=3, include_intercept=False)[4] 3.524e+05 1.97e+04 17.854 0.000 3.14e+05 3.91e+05
bs(train, knots=(3, 6, 9), degree=3, include_intercept=False)[5] 4.232e+05 4.12e+04 10.269 0.000 3.42e+05 5.04e+05
# obtain predicted value from model spline
train_dataset2['pred_spline'] = spline1.predict()

# obtain the residuals from statsmodel (resid)
train_dataset2['error_2_spline'] = train_dataset2['median_house_value'] - train_dataset2['pred_spline']

# square the residuals 
train_dataset2['error_sq_2_spline'] = train_dataset2['error_2_spline']**2

# show df
train_dataset2.head(5)
longitude latitude housing_median_age total_rooms total_bedrooms population households median_income median_house_value ocean_proximity ... bedrooms_per_rooms average error error_sq pred error_2 error_sq_2 pred_spline error_2_spline error_sq_2_spline
14185 -117.08 32.70 37 2176 418 1301 375 2.8750 98900.0 NEAR COAST ... 0.192096 196674.726944 -97774.726944 9.559897e+09 157722.078408 -58822.078408 3.460037e+09 154283.626628 -55383.626628 3.067346e+09
6125 -117.91 34.11 20 3158 684 2396 713 3.5250 153000.0 NEAR COAST ... 0.216593 196674.726944 -43674.726944 1.907482e+09 185009.672924 -32009.672924 1.024619e+09 183299.913013 -30299.913013 9.180847e+08
14095 -117.10 32.75 11 2393 726 1905 711 1.3448 91300.0 NEAR COAST ... 0.303385 196674.726944 -105374.726944 1.110383e+10 119491.512237 -28191.512237 7.947614e+08 96456.816818 -5156.816818 2.659276e+07
14359 -117.22 32.74 52 1260 202 555 209 7.2758 345200.0 NEAR COAST ... 0.160317 196674.726944 148525.273056 2.205976e+10 394457.593271 -49257.593271 2.426310e+09 372578.112664 -27378.112664 7.495611e+08
18004 -121.99 37.29 32 2930 481 1336 481 6.4631 344100.0 NEAR COAST ... 0.164164 196674.726944 147425.273056 2.173421e+10 334928.406586 9171.593414 8.411813e+07 319906.620313 24193.379687 5.853196e+08

5 rows × 23 columns

# Sum of squared residuals (SSR)
SSR_spline = train_dataset2['error_sq_2_spline'].sum()
print('Sum of squared residuals (SSR) of Cubic Spline:',SSR_spline)

# Explained sum of squares  (SS_M = TSS - SS_R)
SSM_spline = TSS - SSR_spline
print('SSM of Cubic Spline: ', SSM_spline)

# R_Squared: explained sum of squared residuals
R_squared_spline = SSM_spline / TSS
print('R squared of Cubic Spline:', R_squared_spline)
Sum of squared residuals (SSR) of Cubic Spline: 71524714796268.75
SSM of Cubic Spline:  99017711465643.12
R squared of Cubic Spline: 0.5806045664764725
print(skm.r2_score(train_dataset2['median_house_value'], train_dataset2['pred_spline']))
0.5806045664764725

R squared ist schlechter geworden als in den vorherigen Modellen. Allerdings liegt dies darin begründet, dass nur ein Feature verwendet wurde. Ein direkter Vergleich ist also nicht möglich.

# Create observations
xp = np.linspace(X_test.min(),X_test.max(), 100)
# Make some predictions
pred = spline1.predict(dmatrix("bs(xp, knots=(3, 6, 9), include_intercept=False)", {"xp": xp}, return_type='dataframe'))

# plot
sns.scatterplot(x=X_train, y=y_train)

plt.plot(xp, pred, label='Cubic spline with degree=3 (3 knots)', color='orange')
plt.legend();
_images/regression_153_0.png

Model 4 - Natural Spline

transformed_x3 = dmatrix("cr(train,df = 3)", {"train": X_train}, return_type='dataframe')

spline2 = sm.GLM(y_train, transformed_x3).fit()
# obtain predicted value from model spline
train_dataset2['pred_spline'] = spline2.predict()

# obtain the residuals from statsmodel (resid)
train_dataset2['error_2_spline'] = train_dataset2['median_house_value'] - train_dataset2['pred_spline']

# square the residuals 
train_dataset2['error_sq_2_spline'] = train_dataset2['error_2_spline']**2

# show df
train_dataset2.head(5)
longitude latitude housing_median_age total_rooms total_bedrooms population households median_income median_house_value ocean_proximity ... bedrooms_per_rooms average error error_sq pred error_2 error_sq_2 pred_spline error_2_spline error_sq_2_spline
14185 -117.08 32.70 37 2176 418 1301 375 2.8750 98900.0 NEAR COAST ... 0.192096 195254.121813 -96354.121813 9.284117e+09 157516.967720 -58616.967720 3.435949e+09 153409.972829 -54509.972829 2.971337e+09
6125 -117.91 34.11 20 3158 684 2396 713 3.5250 153000.0 <1H OCEAN ... 0.216593 195254.121813 -42254.121813 1.785411e+09 182232.226163 -29232.226163 8.545230e+08 180983.759544 -27983.759544 7.830908e+08
14095 -117.10 32.75 11 2393 726 1905 711 1.3448 91300.0 NEAR COAST ... 0.303385 195254.121813 -103954.121813 1.080646e+10 122387.245306 -31087.245306 9.664168e+08 90298.835046 1001.164954 1.002331e+06
14359 -117.22 32.74 52 1260 202 555 209 7.2758 345200.0 NEAR COAST ... 0.160317 195254.121813 149945.878187 2.248377e+10 394045.429765 -48845.429765 2.385876e+09 352306.649672 -7106.649672 5.050447e+07
18004 -121.99 37.29 32 2930 481 1336 481 6.4631 344100.0 <1H OCEAN ... 0.164164 195254.121813 148845.878187 2.215510e+10 331461.374605 12638.625395 1.597349e+08 313580.198527 30519.801473 9.314583e+08

5 rows × 23 columns

# Sum of squared residuals (SSR)
SSR_spline = train_dataset2['error_sq_2_spline'].sum()
print('Sum of squared residuals (SSR) of Natural Spline:',SSR_spline)

# Explained sum of squares  (SS_M = TSS - SS_R)
SSM_spline = TSS - SSR_spline
print('SSM of Natural Spline: ', SSM_spline)

# R_Squared: explained sum of squared residuals
R_squared_spline = SSM_spline / TSS
print('R squared of Natural Spline:', R_squared_spline)
Sum of squared residuals (SSR) of Natural Spline: 67847658465067.86
SSM of Natural Spline:  95453782865582.95
R squared of Natural Spline: 0.5845250482040099

R-squared ist im Vergleich zur Cubic Spline nochmal gesunken. Auch hier wird im Modell nur eine unabhängige Variable verwendet.

# Make predictions
pred = spline2.predict(dmatrix("cr(xp, df=3)", {"xp": xp}, return_type='dataframe'))

# plot
sns.scatterplot(x=X_train, y=y_train)
plt.plot(xp, pred, color='orange', label='Natural spline with df=3')
plt.legend()
<matplotlib.legend.Legend at 0x1a182d14490>
_images/regression_159_1.png

Scikit Learn

Data Preprocessing Pipeline

Initialisierung der Preprocessing Pipeline. Mithilfe einer Preprocessing Pipelines müssen die Daten nicht manuell transformiert werden, da diesse Aufgaben von der Pipeline übernommen werden.
Für numerische Variablen gehört dazu die Normalisierung der Variablen, sowie die Belegung nicht befüllter Datensätze mit dem Durchschnittswert des Features.
Bei kategorialen Variablen werden fehlnde Werte durch eine Konstante belegt und automatisiert Dummy Variablen für die verschiedenen Ausprägungen der Variablen erzeugt.

# for numeric features
numeric_transformer = Pipeline(steps=[
    ('imputer', SimpleImputer(strategy='median')),
    ('scaler', StandardScaler()) #Normalisierung (X=(X-Mittelwert) / Standardabweichung)
    ])
# for categorical features  
categorical_transformer = Pipeline(steps=[
    ('imputer', SimpleImputer(strategy='constant', fill_value='missing')), #Konstante bei fehlenden Werten reinschreiben
    ('onehot', OneHotEncoder(handle_unknown='ignore')) #pro Featureausprägung eine Spalte, zutreffendes Feature hat dann Wert 1 und die anderen Wert 0
    ])
# Pipeline
preprocessor = ColumnTransformer(transformers=[
    ('num', numeric_transformer, selector(dtype_exclude="category")),
    ('cat', categorical_transformer, selector(dtype_include="category"))
        ])

Datenvorbereitung

Für Scikit Learn muss die abhängige Variable von den unabhängigen getrennt werden. Als Datenbasis wird hier auf train_dataset2 zugegriffen. In diesem Datenset wurden die Outlier nach Cook’s Distance bereits entfernt. Als Features sollen die durch die Filter Methode identifizierten und bereits in statsmodels verwendeten Features genutzt werden.

# create label y
y = train_dataset2['median_house_value']

# create features
X= train_dataset2[['median_income', 'ocean_proximity', 'bedrooms_per_rooms', 'people_per_household', 'housing_median_age']]

Scikit Learn bietet eine eigene Methode zum Datensplitting. Das Datenset train_dataset2 aus statsmodels wird nun nochmal in Test- und Trainingsdaten aufgeteilt. Dadurch ist Validierung der Modelle bereits mit diesem Testset möglich, bevor am Ende dieses Notebooks mit den bereits zu Beginn zur Seite gelegten Testdaten das beste Modell evaluiert wird.

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=10)
print("Training size:", len(X_train))
print("Test size:", len(X_test))
Training size: 12404
Test size: 3102

Model 5 - Simple Regression

Nun wird ein Pipeline Object erzeugt, welches zunächst die Data Preprocessing Pipeline ausführt und anschließend das Regressionsmodell.

# Create pipeline with model
lm_pipe = Pipeline(steps=[
    ('preprocessor', preprocessor),
    ('lm', LinearRegression())
                        ])

# Fit model
lm_pipe.fit(X_train, y_train)
Pipeline(steps=[('preprocessor',
                 ColumnTransformer(transformers=[('num',
                                                  Pipeline(steps=[('imputer',
                                                                   SimpleImputer(strategy='median')),
                                                                  ('scaler',
                                                                   StandardScaler())]),
                                                  <sklearn.compose._column_transformer.make_column_selector object at 0x0000019AF0E691C0>),
                                                 ('cat',
                                                  Pipeline(steps=[('imputer',
                                                                   SimpleImputer(fill_value='missing',
                                                                                 strategy='constant')),
                                                                  ('onehot',
                                                                   OneHotEncoder(handle_unknown='ignore'))]),
                                                  <sklearn.compose._column_transformer.make_column_selector object at 0x0000019AF0E69700>)])),
                ('lm', LinearRegression())])

In der Variable features_names werden die Namen der Features gespeichert.

feature_names = np.concatenate((X_train.drop(['ocean_proximity'], axis=1).columns.to_numpy(), lm_pipe.named_steps['preprocessor'].transformers_[1][1]['onehot'].get_feature_names_out()))
feature_names
array(['median_income', 'bedrooms_per_rooms', 'people_per_household',
       'housing_median_age', 'x0_INLAND', 'x0_NEAR COAST'], dtype=object)

Der eben erzeugte feature_names Array kann nun mit den Koeffizienten des Modells ergänzt werden, um die Koeffizienten je Feature auszugeben.

# Obtain model coefficients and names
print('Koeffizienten und Namen: \n', list(zip(lm_pipe.named_steps['lm'].coef_, feature_names)))
Koeffizienten und Namen: 
 [(85229.2089681146, 'median_income'), (19379.143797351342, 'bedrooms_per_rooms'), (-19642.68963026445, 'people_per_household'), (11314.894820139021, 'housing_median_age'), (-28646.670935754682, 'x0_INLAND'), (28646.670935754682, 'x0_NEAR COAST')]

Das trainierte Modell wird nun für die Berechnung der median_house_values der Trainingsdaten verwendet.

y_pred = lm_pipe.predict(X_train)
y_pred
array([186053.02255921,  51495.90700283, 263988.77300622, ...,
       152104.43702161, 122857.34230185, 120051.11769621])
# Metrics trainings set
print('R squared training set:', round(skm.r2_score(y_train, y_pred)*100, 2))
print('MSE training set:', round(skm.mean_squared_error(y_train, y_pred), 2))
print('RMSE training set:', round(skm.mean_squared_error(y_train, y_pred, squared=False), 2))
R squared training set: 74.18
MSE training set: 2816155847.98
RMSE training set: 53067.47

Das Modell hat einen ähnlichen R-squared wie die OLS Regression in statsmodels. Der Root Mean Squared Error des Modells beträgt 50968,88$. Das bedeutet, dass der eigentliche Wert von median_house_value im Schnitt um diesen Betrag von dem durch das Modell berechneten Ergebnis abweicht.

Im Folgenden wird das Modell mit den im Scikit Learn Datensplit erzeugten Validierungsdaten getestet.

# Validation data
pred = lm_pipe.predict(X_test)
r2_test = round(skm.r2_score(y_test, pred)*100, 2)
print('R squared validation set:', r2_test)
mse_test =skm.mean_squared_error(y_test, pred)
print('MSE validation set', round(mse_test, 2))
print('RMSE training set:', round(skm.mean_squared_error(y_test, pred, squared=False), 2))
R squared validation set: 74.67
MSE validation set 2877367258.45
RMSE training set: 53641.1

Der R-squared Wert hat sich mit den Validierungsdaten verbessert. Dies ist sehr ungewöhnliche, da unter Verwendung unbekannter Variablen die Metriken normalerweise etwas schlechtere Werte einnehmen. Dies könnte allerdings darin begründet sein, dass mit diesen Validierungsdaten bereits EDA, Feature Engineering und auch Cook’s Distance Outlier Entfernung durchgeführt wurde.
Der durchschnittliche Fehler MSE und auch die durchschnittliche Abweichung von berechneten zu tatsächlichen Modellen sind etwas schlechter geworden. Allerdings nur geringfügig. Das Modell wird daher als gut erachtet.

# get absolute values of coefficients
importance = np.abs(lm_pipe.named_steps['lm'].coef_)

sns.barplot(x=importance, 
            y=feature_names);
_images/regression_186_0.png

Model 6 - Lasso

vorgegebenes Alpha

In diesem Unterkapitel wird eine Lasso Regression für das Modell durchgeführt. Der Wert für alpha wird initial mit 1 belegt.

# Create pipeline with model
lasso_pipe = Pipeline(steps=[
    ('preprocessor', preprocessor),
    ('lasso', Lasso(alpha=1, max_iter=10000))
                        ])

lasso_pipe.fit(X_train, y_train)
Pipeline(steps=[('preprocessor',
                 ColumnTransformer(transformers=[('num',
                                                  Pipeline(steps=[('imputer',
                                                                   SimpleImputer(strategy='median')),
                                                                  ('scaler',
                                                                   StandardScaler())]),
                                                  <sklearn.compose._column_transformer.make_column_selector object at 0x0000019AF0E691C0>),
                                                 ('cat',
                                                  Pipeline(steps=[('imputer',
                                                                   SimpleImputer(fill_value='missing',
                                                                                 strategy='constant')),
                                                                  ('onehot',
                                                                   OneHotEncoder(handle_unknown='ignore'))]),
                                                  <sklearn.compose._column_transformer.make_column_selector object at 0x0000019AF0E69700>)])),
                ('lasso', Lasso(alpha=1, max_iter=10000))])
# Metrics trainings set
pred_train = lasso_pipe.predict(X_train)
print('R squared training set:', round(skm.r2_score(y_train, pred_train)*100, 2))
print('MSE training set:', round(skm.mean_squared_error(y_train, pred_train), 2))
print('RMSE training set:', round(skm.mean_squared_error(y_train, pred_train, squared=False), 2))
R squared training set: 74.18
MSE training set: 2816155856.76
RMSE training set: 53067.47

R-squared der Lasso Regression hat sich im Vergleich zur Simple Regression im Modell zuvor nicht verändert. Auch MSE und RMSE sind gleich geblieben. Dies ist auf den initial gewählten, sehr kleinen Wert für alpha von 1 zurückzuführen. Ein geringer Wert in Alpha verändert die Höhe der Koeffizienten im Modell kaum. Daher bleiben auch die Berechnungen des Modells unverändert und somit die Metriken.

# Metrics validation set
pred = lasso_pipe.predict(X_test)
print('R squared validation set:', round(skm.r2_score(y_test, pred)*100, 2))
print('MSE validation set:', round(skm.mean_squared_error(y_test, pred), 2))
print('RMSE validation set:', round(skm.mean_squared_error(y_test, pred, squared=False), 2))
R squared validation set: 74.67
MSE validation set: 2877369775.47
RMSE validation set: 53641.12

Auch die Metriken des Validation Set haben sich nicht verändert im Vergleich zur Simple Regression.

# Obtain model coefficients and names
print('Koeffizienten und Namen: \n', list(zip(lasso_pipe.named_steps['lasso'].coef_, feature_names)))
Koeffizienten und Namen: 
 [(85226.87228779327, 'median_income'), (19376.806990272275, 'bedrooms_per_rooms'), (-19641.750378420238, 'people_per_household'), (11314.128700828936, 'housing_median_age'), (-57290.99849238739, 'x0_INLAND'), (2.77139978999994e-10, 'x0_NEAR COAST')]
# get absolute values of coefficients
importance = np.abs(lasso_pipe.named_steps['lasso'].coef_)

sns.barplot(x=importance, 
            y=feature_names);
_images/regression_196_0.png

Das Diagramm visualiert die Feature Importance. Je länger der Balken, desto eine höhere Gewichtung bekommt das Feature im Modell. Die Feature Importance ergibt sich auch den Koeffizienten der Features. median_income ist das wichtigste Feature im Modell, wie an dem großen Balken sowie dem Koeffizienten des Features zu erkennen.
Die Ausprägung NEAR COAST der Variablen ocean_proximity wurde von der Lasso Regression aus der Regression entfernt. Dieses Verhalten ist bei Lasso gewünscht. Durch die L1-Regularisierung mit alpha führt der Algorithmus automatisch eine Feature Selection durch. Dabei werden relevante Features durch höhere Koeffizienten stärker gewichtet und unwichtige Features durch kleinere Koeffizienten schäwcher gewichtet bzw. sogar komplett aus dem Modell entfernt, indem der Koeffizient auf den Wert 0 reduziert wird. Durch diese Methode soll ein overfitting des Modells durch zu viele Variablen verhindert werden.

k-fold cross validation zur Ermittlung von alpha

from sklearn.linear_model import LassoCV

# Lasso with 5 fold cross-validation
model_pipe = Pipeline(steps=[
    ('preprocessor', preprocessor),
    ('lassoCV', LassoCV(cv=5, random_state=0, max_iter=10000))
                        ])
# Fit model
model_pipe.fit(X_train, y_train)
Pipeline(steps=[('preprocessor',
                 ColumnTransformer(transformers=[('num',
                                                  Pipeline(steps=[('imputer',
                                                                   SimpleImputer(strategy='median')),
                                                                  ('scaler',
                                                                   StandardScaler())]),
                                                  <sklearn.compose._column_transformer.make_column_selector object at 0x0000019AF0E691C0>),
                                                 ('cat',
                                                  Pipeline(steps=[('imputer',
                                                                   SimpleImputer(fill_value='missing',
                                                                                 strategy='constant')),
                                                                  ('onehot',
                                                                   OneHotEncoder(handle_unknown='ignore'))]),
                                                  <sklearn.compose._column_transformer.make_column_selector object at 0x0000019AF0E69700>)])),
                ('lassoCV', LassoCV(cv=5, max_iter=10000, random_state=0))])
#Show best value of penalization chosen by cross validation:
alpha = model_pipe.named_steps['lassoCV'].alpha_
print('Best alpha chosen by cross validation:', alpha)
Best alpha chosen by cross validation: 78.94891231105018

Auf Basis von k-fold cross validation wurde der beste Wert für Alpha ermittelt. Dieser liegt bei 79,59. Bei k-fold cross validation wird das Trainingsdatenset nochmal in k kleinere Datensets unterteilt. Eines Dieser kleineren Datensets wird zur Seite gelegt und das Modell mithilfe der anderen und einem Wert für Alpha trainiert. Anschließend wird mit dem zur Seite gelegten kleinen Datenset eine Validierung des Modells mit den unbekannten Daten durchgeführt.
Dieser Vorgang wird für jedes der k - Datensets für verschiedene Werte von Alpha wiederholt. Das Alpha mit dem besten Ergebnis bei der Validierung wird hier ausgegeben.

Mit dem zuvor ermittelten, besten Wert von Alpha wird nun noch ein finales Lasso Modell auf Basis aller Trainingsdaten trainiert.

# Set best alpha
lasso_best_pipe = Pipeline(steps=[
    ('preprocessor', preprocessor),
    ('lasso_best', Lasso(alpha=alpha))
                        ])
lasso_best_pipe.fit(X_train, y_train)
lasso_best = lasso_best_pipe.named_steps['lasso_best']
print('Model coefficients and names: \n', list(zip(lasso_best.coef_, feature_names)))
Model coefficients and names: 
 [(85045.084637978, 'median_income'), (19194.892651470986, 'bedrooms_per_rooms'), (-19568.52049548915, 'people_per_household'), (11254.486632101924, 'housing_median_age'), (-57108.04741004523, 'x0_INLAND'), (1.98856656715685e-10, 'x0_NEAR COAST')]

Durch den größeren Wert von Alpha wurden die Koeffizienten der einzelnen Features etwas verkleinert. median_income hatte zuvor beispielsweise einen Wert von 85260,00. Nun liegt der Wert bei 85079,34.

# get absolute values of coefficients
importance = np.abs(lasso_best.coef_)

sns.barplot(x=importance, 
            y=feature_names);
_images/regression_206_0.png

An der Reihenfolge der Feature Importance hat sich durch die Veränderung von Alpha wenig verändert.

# Metrics trainings set
pred_train = lasso_best_pipe.predict(X_train)
print('R squared training set:', round(skm.r2_score(y_train, pred_train)*100, 2))
print('MSE training set:', round(skm.mean_squared_error(y_train, pred_train), 2))
print('RMSE training set:', round(skm.mean_squared_error(y_train, pred_train, squared=False), 2))
R squared training set: 74.18
MSE training set: 2816210527.95
RMSE training set: 53067.98
# Metrics tests set
pred = lasso_best_pipe.predict(X_test)
print('R squared validation set:', round(skm.r2_score(y_test, pred)*100, 2))
print('MSE validation set:', round(skm.mean_squared_error(y_test, pred), 2))
print('RMSE validation set:', round(skm.mean_squared_error(y_test, pred, squared=False), 2))
R squared validation set: 74.67
MSE validation set: 2877619408.17
RMSE validation set: 53643.45

Da sich die Koeffizienten durch den neuen Wert von Alpha kaum verändert haben, sind auch die Kennzahlen zur Performance Messung des Modells unverändert.

Model 7 - Natural Splines

In diesem Modell wird eine scikit learn spline in Kombination mit einer Data Preprocessing Pipeline trainiert.

# use a spline wit 4 knots and 3 degrees with a ridge regressions
spline = make_pipeline(SplineTransformer(n_knots=4, degree=3), 
                       Ridge(alpha=1))

# Integrate Preprocessor Pipeline
spline_pipe = Pipeline(steps=[
    ('preprocessor', preprocessor),
    ('spline', spline)
                        ])
                     
spline_pipe.fit(X_train, y_train)
Pipeline(steps=[('preprocessor',
                 ColumnTransformer(transformers=[('num',
                                                  Pipeline(steps=[('imputer',
                                                                   SimpleImputer(strategy='median')),
                                                                  ('scaler',
                                                                   StandardScaler())]),
                                                  <sklearn.compose._column_transformer.make_column_selector object at 0x0000019AF0E691C0>),
                                                 ('cat',
                                                  Pipeline(steps=[('imputer',
                                                                   SimpleImputer(fill_value='missing',
                                                                                 strategy='constant')),
                                                                  ('onehot',
                                                                   OneHotEncoder(handle_unknown='ignore'))]),
                                                  <sklearn.compose._column_transformer.make_column_selector object at 0x0000019AF0E69700>)])),
                ('spline',
                 Pipeline(steps=[('splinetransformer',
                                  SplineTransformer(n_knots=4)),
                                 ('ridge', Ridge(alpha=1))]))])
# Metrics trainings set
pred_train = spline_pipe.predict(X_train)
print('R squared training set:', round(skm.r2_score(y_train, pred_train)*100, 2))
print('MSE training set:', round(skm.mean_squared_error(y_train, pred_train), 2))
print('RMSE training set:', round(skm.mean_squared_error(y_train, pred_train, squared=False), 2))
R squared training set: 75.09
MSE training set: 2716528957.34
RMSE training set: 52120.33
# Metrics tests set
pred = spline_pipe.predict(X_test)
print('R squared validation set:', round(skm.r2_score(y_test, pred)*100, 2))
print('MSE validation set:', round(skm.mean_squared_error(y_test, pred), 2))
print('RMSE validation set:', round(skm.mean_squared_error(y_test, pred, squared=False), 2))
R squared validation set: 75.52
MSE validation set: 2781732015.02
RMSE validation set: 52742.13

Durch die Verwendung der Spline hat sich das Modell in Bezug auf R-sqared und auch RMSE nochmal verbessert. Das Modell deckt nun knapp 76% der Varianz ab und die durchschnittliche Abweichung der tatsächlichen Werten zu den berechneten Werten beträgt knapp 51000$.

# Create observations
xp = pd.DataFrame(data=np.linspace(X_test.min(),X_test.max(), 100), columns = X_train.drop(['ocean_proximity'], axis=1).columns)
xp = xp.assign(ocean_proximity=lambda xp: 'INLAND')


# Make some predictions
pred_new = spline_pipe.predict(xp)

# plot
sns.scatterplot(x=X_train['median_income'], y=y_train)

plt.plot(xp['median_income'], pred_new, label='Cubic spline with degree=3', color='orange')
plt.legend();
C:\Users\Admin\AppData\Local\Temp/ipykernel_9496/3121680389.py:2: FutureWarning:

Dropping of nuisance columns in DataFrame reductions (with 'numeric_only=None') is deprecated; in a future version this will raise TypeError.  Select only valid columns before calling the reduction.
_images/regression_217_1.png

Evaluierung der besten 2 Modelle mit Test Datenset

In diesem Kapitel wird eine Evaluierung der zwei besten Modellen mithilfe der zu Beginn beiseite gelegten Testdaten durchgeführt.
Zunächst müssen die Testdaten transformiert werden, sodass sie in gleicher Form wie die Trainingsdaten vorliegen und von den Modellen verwendet werden können.

# write out a dict with the mapping of old to new
remap_cat_dict = {
    'NEAR OCEAN': 'NEAR COAST',
    'NEAR BAY': 'NEAR COAST',
    'INLAND': 'INLAND',
    '<1H OCEAN': 'NEAR COAST',
    'ISLAND' : 'NEAR COAST'}

test_dataset.ocean_proximity = test_dataset.ocean_proximity.map(remap_cat_dict).astype('category')
y_test = test_dataset['median_house_value']

X_test = test_dataset[['median_income', 'ocean_proximity', 'bedrooms_per_rooms', 'people_per_household', 'housing_median_age']]

Evaluation Statsmodels OLS lm2

# Metrics tests set
X_test = X_test.assign(pred_test1=lambda X_test: lm2.predict(X_test))
r2_test1 = round(skm.r2_score(y_test, X_test.pred_test1)*100, 2)
print('R squared Statsmodels OLS with Test Set:', r2_test1)
print('MSE Statsmodels OLS with Test Set:', round(skm.mean_squared_error(y_test, X_test.pred_test1), 2))
print('RMSE Statsmodels OLS with Test Set:', round(skm.mean_squared_error(y_test, X_test.pred_test1, squared=False), 2))
R squared Statsmodels OLS with Test Set: -2.13
MSE Statsmodels OLS with Test Set: 13381618868.21
RMSE Statsmodels OLS with Test Set: 115678.95

R-squared ist in Model lm2 negativ. Somit generalisiert das Modell schlecht auf unbekannte Testdaten und stellt sogar eine Verschlechterung im Vergleich zur Verwendung des Durchschnitts als Modell dar.
Grund dafür könnte ein Overfitting während dem Modellierungsprozess sein.
Root Mean Squared Error und somit eine durchschnittliche Abweichung von tatsächlichem und berechnetem median_house_value in Höhe von 115678,95$ ist ebenso kein zufriedenstellendes Ergebnis.

Die berechneten Werte sollen in einem Scatterplot den tatsächlichen Werten gegenübergestellt werden. Aufgrund der Datenmenge sollen nur einige Datensätze visualisiert werden. Mithilfe eines Data Splittings werden zufällig einige Datenpunkte ausgewählt.

X_notplot, X_plot, y_notplot, y_plot = train_test_split(X_test, y_test, test_size=0.1, random_state=10)
# plot
sns.scatterplot(x=X_plot['median_income'], y=y_plot, label='actual')

plt.plot(X_plot['median_income'], X_plot['pred_test1'], label='predicted', color='orange', marker='.', linestyle='')
plt.legend();
_images/regression_228_0.png

Im Plot ist zu erkennen, dass bei bei einem median_income von 2 einen Datenpunkt gibt, dessen median_house_value negativ geschätzt wurde. Dieser Wert ist unplausibel, da ein Haus keinen negativen Wert hat.
Den Koeffizienten des Modells lm2 ist zu entnehmen, dass nur das Feature people_per_household eine negative Abhängigkeit auf median_house_value hat. Die Ursache für einen negativen Wert wird also durch einen hohen Wert im Feature people_per_household verursacht.

lm2.params
Intercept                        -55841.725014
ocean_proximity[T.NEAR COAST]     57312.739989
median_income                     48796.057520
bedrooms_per_rooms               368968.477487
people_per_household             -25917.607221
housing_median_age                  924.724553
dtype: float64
test_dataset.loc[test_dataset['people_per_household'] > 100]
longitude latitude housing_median_age total_rooms total_bedrooms population households median_income median_house_value ocean_proximity people_per_household rooms_per_household bedrooms_per_people bedrooms_per_rooms
13034 -121.15 38.69 52 240 44 6675 29 6.1359 225000.0 INLAND 230.172414 8.275862 0.006592 0.183333

Tatsächlich existiert ein Datensatz mit über 100 people_per_household. Dieser ist vermutlich für den negativen R-squared Wert verantwortlich. Daher wird dieser testweise entfernt und R-squared erneut berechnet.

X_test2 = X_test.drop(13034)
y_test2 = y_test.drop(13034)
# Metrics tests set
X_test2 = X_test2.assign(pred_test1=lambda X_test2: lm2.predict(X_test2))
r2_test1 = round(skm.r2_score(y_test2, X_test2.pred_test1)*100, 2)
print('R squared Statsmodels OLS with Test Set:', r2_test1)
print('MSE Statsmodels OLS with Test Set:', round(skm.mean_squared_error(y_test2, X_test2.pred_test1), 2))
print('RMSE Statsmodels OLS with Test Set:', round(skm.mean_squared_error(y_test2, X_test2.pred_test1, squared=False), 2))
R squared Statsmodels OLS with Test Set: 61.37
MSE Statsmodels OLS with Test Set: 5063025615.54
RMSE Statsmodels OLS with Test Set: 71154.94

Wie vermutet, lag der Grund für den negativen Wert in diesem Ausreiser. Jedoch darf auf den Testdaten gemäß der Literatur keine EDA oder Feature Engineering betrieben werden.
Das Modell lm2 ist daher kein gutes Modell, welches ausreichend unbekannte Daten generalisiert.

Evaluation Scikit Learn Spline

#Calculate R2
X_test = X_test.assign(pred_test2=lambda X_test: spline_pipe.predict(X_test))
r2_test = round(skm.r2_score(y_test, X_test.pred_test2)*100, 2)
print('R squared test set:', r2_test)
mse_test =skm.mean_squared_error(y_test, X_test.pred_test2)
print('MSE test set', round(mse_test, 2))
print('RMSE Statsmodels OLS with Test Set:', round(skm.mean_squared_error(y_test, X_test.pred_test2, squared=False), 2))
R squared test set: 65.25
MSE test set 4552968088.14
RMSE Statsmodels OLS with Test Set: 67475.69

Dieses Modell berechnet aus dem kompletten Test Datenset einen zufriedenstellenden Wert für R-squared.
Auch RMSE mit einem Wert von 67475,69$ ist in Ordnung.

# plot
sns.scatterplot(x=X_plot['median_income'], y=y_plot, label='actual')

plt.plot(X_plot['median_income'], X_plot['pred_test2'], label='predicted', color='orange', marker='.', linestyle='')
plt.legend();
_images/regression_239_0.png

Den Analysen entsprechend ist das beste Modell ist die Spline von Scikit Learn mit einem R-squared der Testdaten von 65,25%.
Gemäß dem Data-Science Lifecycle würde dieses Modell jetzt deployed werden. Dieser Schritt wird allerdings nicht mehr im Rahmen dieses Projekts durchgeführt.